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Abstract 

DDSCAT 6.1 is a freely available software package which applies the "discrete dipole approxima- 
tion" (DDA) to calculate scattering and absorption of electromagnetic waves by targets with arbitrary 
geometries and complex refractive index. The DDA approximates the target by an array of polarizable 
points. DDSCAT 6.1 allows accurate calculations of electromagnetic scattering from targets with "size 
parameters" 27racff/ A < 15 provided the refractive index m is not large compared to unity (|to— 1| < 2). 
DDSCAT 6.1 includes the option of using the FFTW (Fastest Fourier Transform in the West) package. 
DDSCAT 6.1 also includes support for MP I (Message Passing Interface), permitting parallel calculations 
on multiprocessor systems. 

We also make available a "plain" distribution of DDSCAT 6.1 that does not include support for MP I, 
FFTW, or netCDF, but is much simpler to install than the full distribution. 

The DDSCAT package is written in Fortran and is highly portable. The program supports calcu- 
lations for a variety of target geometries (e.g., ellipsoids, regular tetrahedra, rectangular solids, finite 
cyhnders, hexagonal prisms, etc.). Target materials may be both inhomogeneous and anisotropic. It is 
straightforward for the user to "import" arbitrary target geometries into the code, and relatively straight- 
forward to add new target generation capability to the package. DDSCAT automatically calculates total 
cross sections for absorption and scattering and selected elements of the Mueller scattering intensity 
matrix for specified orientation of the target relative to the incident wave, and for specified scattering 
directions. 

This User Guide explains how to use DDSCAT 6.1 to carry out electromagnetic scattering calcula- 
tions. CPU and memory requirements are described. 
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1 Introduction 

DDSCAT 6.1 is a Fortran software package to calculate scattering and absorption of electromagnetic 
waves by targets with arbitrary geometries using the "discrete dipole approximation" (DDA). In this 
approximation the target is replaced by an array of point dipoles (or, more precisely, polarizable points); 
the electromagnetic scattering problem for an incident periodic wave interacting with this array of point 
dipoles is then solved essentially exactly. The DDA (sometimes referred to as the "coupled dipole 
approximation") was apparently first proposed by Purcell & Pennypacker (1973). DDA theory was 
reviewed and developed further by Draine (1988), Draine & Goodman (1993), and recently reviewed by 
Draine & Flatau (1994) and Draine (2000). 

DDSCAT 6.1 is a Fortran implementation of the DDA developed by the authorsQ 
It is intended to be a versatile tool, suitable for a wide variety of applications ranging from interstellar 
dust to atmospheric aerosols. As provided, DDSCAT 6.1 should be usable for many applications without 
modification, but the program is written in a modular form, so that modifications, if required, should be 
fairly straightforward. 

The authors make this code openly available to others, in the hope that it will prove a useful tool. We 
ask only that; 

• If you publish results obtained using DDSCAT, please acknowledge the source of the code. 

• If you discover any errors in the code or documentation, please promptly communicate them to the 
authors. 

• You comply with the "copyleft" agreement (more formally, the GNU General Public License) of 
the Free Software Foundation: you may copy, distribute, and/or modify the software identified as 
coming under this agreement. If you distribute copies of this software, you must give the recip- 
ients all the rights which you have. See the file doc/copyleft distributed with the DDSCAT 
software. 

We also strongly encourage you to send email to the authors identifying yourself as a user of DDSCAT; 
this will enable the authors to notify you of any bugs, corrections, or improvements in DDSCAT. 

The current version, DDSCAT 6.1, uses the DDA formulae from Draine (1988), with dipole po- 
larizabilities determined from the Lattice Dispersion Relation (Draine & Goodman 1993). The code 
incorporates Fast Fourier Transform (FFT) methods (Goodman, Draine, & Flatau 1991). 

We refer you to the list of references at the end of this document for discussions of the theory and 
accuracy of the DDA [first see the recent reviews by Draine and Flatau (1994) and Draine (2000)]. In 
yHwe describe the principal changes between DDSCAT 6.1 and the previous releases. The succeeding 
sections contain instructions for: 

• compiling and linking the code; 

• running a sample calculation; 

• understanding the output from the sample calculation; 

• modifying the parameter file to do your desired calculations; 

• specifying target orientation; 

' The release history of DDSCAT is as follows: 

• DDSCAT 4b: Released 1993 March 12 

• DDSCAT 4bl: Released 1993 July 9 

• DDSCAT 4c: Although never announced, DDSCAT . 4 c was made available to a number of interested users beginning 1994 
December 18 

• DDSCAT 5a7: Released 1996 

• DDSCAT 5a8: Released 1997 April 24 

• DDSCAT 5a9: Released 1998 December 15 

• DDSCAT SalO: Released 2000 June 15 

• DDSCAT 6.0: Released 2003 September 2 

• DDSCAT 6.1: Released 2004 September 10 
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• changing the DIMENSIONing of the source code to accommodate your desired calculations. 

The instructions for compiling, linking, and running will be appropriate for a UNIX system; slight 
changes will be necessary for non-UNIX sites, but they are quite minor and should present no difficulty. 

Finally, the current version of this User Guide can be obtained from 
|http : / / arxiv . org/abs/astro-ph / 0008151| - you will be offered the options of downloading 
either Postscript or PDF versions of the User Guide. 



2 Applicability of the DDA 

The principal advantage of the DDA is that it is completely flexible regarding the geometry of the target, 
being limited only by the need to use an interdipole separation d small compared to (1) any structural 
lengths in the target, and (2) the wavelength A. Numerical studies (Draine & Goodman 1993; Draine & 
Flatau 1994; Draine 2000) indicate that the second criterion is adequately satisfied if 

|m|fcd<l , (1) 

where m is the complex refractive index of the target material, and k = 2it/X, where A is the wavelength 
in vacuo. However, if accurate calculations of the scattering phase function (e.g., radar or lidar cross 
sections) are desired, a more conservative criterion 

|m|fcd<0.5 (2) 

will ensure that differential scattering cross sections dCsca/dft are accurate to within a few percent of 
the average differential scattering cross section Csca/47r (see Draine 2000). 

Let V be the target volume. If the target is represented by an array of N dipoles, located on a cubic 
lattice with lattice spacing d, then 

V = Nd^ . (3) 
We characterize the size of the target by the "effective radius" 

aeff EE (31//47r)i/3 , (4) 

the radius of an equal volume sphere. A given scattering problem is then characterized by the dimen- 
sionless "size parameter" 

_ 27racff 
X = kacff = — r — . (5) 



A 



The size parameter can be related to N and \m\kd: 



27racff 62.04 / N ^^^^ 



A \m\ \10^ 

Equivalently, the target size can be written 

A /7V^^/^ 



\kd . (6) 



m 



106 



Ocff = 9.873-- — r • \m\kd . (7) 



Practical considerations of CPU speed and computer memory currently available on scientific worksta- 
tions typically limit the number of dipoles employed to < 10® (see iJT6]for limitations on A^ due to 
available RAM); for a given A^, the limitations on \m\kd translate into limitations on the ratio of target 
size to wavelength. 

For calculations of total cross sections Cabs and Cgca, we require \m\kd < 1: 

A / A^ 62.04 / A^ ^ 



m 



Ooff < 9.88-- — or X < — — 77^ • (8) 



106 / \m\ \ 106 
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1 2 3 4 5 6 7 8 9 10 11 12 13 
x = ka 



Fi gure ll Scattering and absorption for a spliere with m = 1.33 + O.Oli. The upper panel shows the exact values of Qsca and 
Qabs, obtained with Mie theory, as functions of a; = ka. The middle and lower panels show fractional errors in Qsca and Qabs, 
obtained using DDSCAT with polarizabilities obtained from the Lattice Dispersion Relation, and labelled by the number N of 
dipoles in each pseudosphere. After Fig. 1 of Draine & Flatau (1994). 



For scattering phase function calculations, we require |m|fc(i < 0.5: 

. A f N 31.02 f nV^^ 

Ooff < 4.94- — - — - or a; < - — - — - . (9) 

\m\ \10^ J \m\ VlO^ / 

It is therefore clear that the DDA is not suitable for very large values of the size parameter x, or very 
large values of the refractive index m. The primary utility of the DDA is for scattering by dielectric 
targets with sizes comparable to the wavelength. As discussed by Draine & Goodman (1993), Draine 
& Flatau (1994), and Draine (2000), total cross sections calculated with the DDA are accurate to a few 
percent provided N > 10'' dipoles are used, criterion ([T]) is satisfied, and the refractive index is not too 
large. 

For fixed \m\kd, the accuracy of the approximation degrades with increasing |m — 1|, for reasons 
having to do with the surface polarization of the target as discussed by Collinge & Draine (2004). With 
the present code, good accuracy can be achieved for |m — 1 1 < 2. 

Examples illustrating the accuracy of the DDA are shown in Figs.[T]-(2]which show overall scattering 
and absorption efficiencies as a function of wavelength for different discrete dipole approximations to a 
sphere, with N ranging from 304 to 59728. The DDA calculations assumed radiation incident along the 
(1,1,1) direction in the "target frame". Figs.[3}S]show the scattering properties calculated with the DDA 
for X ~ ka = 7. Additional examples can be found in Draine & Flatau (1994) and Draine (2000). 



3 DDSCAT 6.1 

3.1 What Does It Calculate? 

DDSCAT 6.1, like previous versions of DDSCAT, solves the problem of scattering and absorption by an 
array of polarizable point dipoles interacting with a monochromatic plane wave incident from infinity. 
DDSCAT 6.1 has the capability of automatically generating dipole array representations for a variety of 
target geometries (see |JT9]) and can also accept dipole array representations of targets supplied by the 
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Figure 3: Differential scattering cross section for m — 1.33 + O.Oli pseudosphere and ka = 7. Lower panel shows fractional 
error compared to exact Mie theory result. The A'^ = 17904 pseudosphere has \m\kd = 0.57, and an rms fractional error in da/dQ 
of 2.4%. After Fig. 5 of Draine & Flatau (1994). 
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Figure 4: Same as Fig.[3]but for m = 2 + i. The A'^ = 59728 pseudosphere has \m\kd = 0.65, and an rms fractional error in 
da/dQ of 6.7%. After Fig. 8 of Draine & Flatau (1994). 



user (although the dipoles must be located on a cubic lattice). The incident plane wave can have arbitrary 
elliptical polarization (see fJlU' ™d the target can be arbitrarily oriented relative to the incident radiation 
(see j fTTI ). The following quantities are calculated by DDSCAT 6.1 : 

• Absorption efficiency factor Qabs = Cabs/'i'aoff' where Cabs is the absorption cross section; 

• Scattering efficiency factor Qsca = C'sca/'ra^g, where Csca is the scattering cross section; 

• Extinction efficiency factor Qoxt = Qsca + Qabs; 

• Phase lag efficiency factor Qpha, defined so that the phase-lag (in radians) of a plane wave after 
propagating a distance L is just TitQphn'^CL^ffL, where rit is the number density of targets. 

• The 4x4 Mueller scattering intensity matrix Sij describing the complete scattering properties of 
the target for scattering directions specified by the user 

• Radiation force efficiency vector Q^ad (see 03. 

• Radiation torque efficiency vector Qr (see ifTS] ). 



3.2 Application to Targets in Dielectric Media 

Let uj be the angular frequency of the incident radiation. For many applications, the target is essen- 
tially in vacuo, in which case the dielectric function e which the user should supply to DDSCAT is the 
actual complex dielectric function Ctargct (i^), or complex refractive index JTitargct (^^) = Retarget of the 
target material. 

However, for many applications of interest (e.g., marine optics, or biological optics) the "target" body 
is embedded in a (nonabsorbing) dielectric medium, with (real) dielectric function emcdium('jj), or (real) 
refractive index mmcdiumC'jj) = Vcmcdium- DDSCAT is fully applicable to these scattering problems, 
except that: 

• The "dielectric function" or "refractive index" supplied to DDSCAT should be the relative dielec- 
tric function 



'^medium V 
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or relative refractive index: 



m{cu) = ^^^^^^^^. (11) 



• The wavelength A specified in ddscat . par should be the wavelength in the medium: 

\= , (12) 

^medium 

where Xuac ~ l-Kcjuj is the wavelength in vacuo. 

The absorption, scattering, extinction, and phase lag efficiency factors Qabs, Qsca, and Qoxt calculated 
by DDSCAT will then be equal to the physical cross sections for absorption, scattering, and extinction 
divided by Tra^jj . For example, the attenuation coefficient for radiation propagating through a medium 
with a number density nt of scatterers will be just a = ntQcxtT^o-cS- Similarly, the phase lag (in radians) 
after propagating a distance L will be ntQphaT^dcS- 

The elements Sij of the 4x4 Mueller scattering matrix S calculated by DDSCAT will be correct for 
scattering in the medium: 



A'' 

2'Kr 



S-I„, (13) 



where Ii„ and Igca are the Stokes vectors for the incident and scattered light (in the medium), r is the 
distance from the target, and A is the wavelength in the medium (eq.[T2b. See |23]for a detailed discussion 
of the Mueller scattering matrix. 

The time-averaged radiative force and torque (see ifTSl l on a target in a dielectric medium are 

Fiad = QprTTfl^ffUiad , (14) 

Trad = QrTraoff^radTr" ' (15) 
ZTT 



where the time-averaged energy density is 



^rad — ^medium ; (16) 



where Eq cos{ujt + cj)) is the electric field of the incident plane wave in the medium. 



4 What's New? 

DDSCAT 6.1 differs from DDSCAT 5a in the following ways: 

1. Via the parameter file ddscat .par, the user can now select up to 9 elements of the Muller 
scattering matrix Sij to be printed out. 

2. The maximum number of iterations allowed has been increased from 300 to 10000, since DDSCAT 
is now increasingly employed for computations that converge relatively slowly: large numbers of 
dipoles, large values of the scattering parameter, or large values of the refractive index. The number 
of iterations used for a solution is recorded in the vivaarbbkccc . sea output files. 

3. In the form distributed, DDSCAT writes only buffered output. When run on a single cpu (with 
MPI disabled), DDSCAT writes "running output" to a file ddscat . log_000. When MPI is 
enabled, cpu nnn writes to a file ddscat . log_nnn, where n=000, 001, 002, etc. 

4. A new target option is available: LYRSLB (a rectangular slab with up to 4 different material layers). 

5. User must now specify (in ddscat . par) which IWAV, IRAD, IWAV to start with. The "default" 
choice will be 0, but when computations have been interupted by, e.g., power outage, there 
will be occasions when it is convenient to be able to resume with the first incomplete calculation. 
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6. A new FFT option is supported: FFTW2 1. This invokes the FFTW ("Fastest Fourier Transform in 
the West") package of Frigo and Johnson. The calls are compatible with versions 2.1.x of FFTW. 
Instructions on compiling and linking to include FFTW are given below ( %.2.5i . 

7. FFT options BRENNR and TMPRTN have been eliminated as they seem to offer no advantages 
relative to GPFAFT or FFTW21. 

8. MP I is now supported. Users can now carry out parallel calculations using DDSCAT 6.1 on 
multiprocessor systems in conformance with the MP I (Message Passing Interface) standards ( i|24] i. 
This of course requires that MP I be akeady installed on the system. 

9. New with DDSCAT 6.1: calculation of (cos^ 6) for the scattered radiation field, in addition to the 
first moment (cos 9) . The format of the output files wwaarbb'kccc .sea has been altered. 

10. New with DDSCAT 6.1: The output files written when the NETCDF option is specified have been 
modified to include more relevant data and to be more easily read. 

11. New with DDSCAT 6.1: The procedure used for averaging over scattering angles has been im- 
proved, both to improve efficiency and to make it easier for the user to manage the balance be- 
tween accuracy and computational burden. Via the parameter file ddscat . par, the user now 
specifies a parameter ETASCA (see ^22]for discussion of the relationship between ETASCA and 
accuracy). Important note: The parameters ICTHM and IPHIM have been ehminated! (See fj8] 
for discussion of the parameter file ddscat . par, or AppendixlAlfor an annotated example of the 
ddscat . par file), ddscat . par files used with previous releases of DDSCAT will not work 
until edited to remove icthm and iphim and add specification of etasca! 

12. New with DDSCAT 6.1: In addition to the full distiibution of DDSCAT 6.1, we also make avail- 
able a "plain" distribution which does not include support for MP I, FFTW, or net CDF, but is much 
simpler for a user to compile and install. This "plain" distribution is recommended for first-time 
users. 

13. New with DDSCAT 6.1: A new target option AN I RFC is available for homogenous rectangular 
targets composed of anisotropic material (see i?T9] l. 

14. New with DDSCAT 6.1: A new target option CYLCAP is available: the target is a cylinder with 
hemispherical end-caps. 

15. 2006.11.28 bug fix: Corrected error in evaluation of P = degree of linear polarization of scattered 
light. Mueller matrix elements Sij were calculated and reported correctly, but P was not correctly 
evaluated. 

5 Obtaining and Installing the Source Code 

5.1 For first-time users: the "plain" distribution 

The "plain" distribution lacks support for MPT, FFTW, netCDF, and reporting of cpu time used in 
different phases of the calculation, but has complete capabilities for target generation and calculating 
scattering using a single CPU. The "plain" distribution is recommended for first-time users. 

5.2 The "full" distribution 

If you require support for MP I, FFTW, or netCDF, you should download the "full" distribution of 
DDSCAT 6.1. 

5.3 Downloading the Source Code 

There are 3 ways to obtain the source code: 

• Point your browser to the DDSCAT 6.1 home page: 
[http : / /www .astro . princeton . edu/^draine/DDSCAT . html| 
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where you will find links for ddscat 6 . l_plain . tgz and ddscat 6 . l_f ull . tgz, gzipped 
tarfiles containg the complete source code for the "plain" and "full" distributions, respectively. 
Right-click on the appropriate link. 

• Point your browser to 

[ftp : / / ftp . astro .princeton . edu/ draine/ scat/ddscat/ver6 . 1| 
right-click on ddscat 6 . l_plain . tgz or ddscat 6 . l_f ull . tgz . 

• Use anonymous ftp to ftp.astro.princeton.edu. 
In response to "Name:" enter anonymous. 

In response to "Password:" enter your full email address. 

cd draine/scat/ddscat/ver6 . 1 
binary ( to enable transfer of binary data 

get ddscat 6 . l_plain . tgz (To get the " plain" distribution) 
or 

get ddscat 6 . l_f ull . tgz (To get the "full" distribution). 

After downloading either ddscat 6 . l_plain . tgz or ddscat 6 . l_f ull . tgz into the directory 
where you would like DDSCAT 6.1 to reside (you should have at least 10 Mbytes of disk space avail- 
able), the source code can be installed as follows: 

5.4 Installing the source code on a Linux/Unix system 

If you are on a Linux system, you should be able to type 

tar xvzf ddscat 6 . l_plain . tgz 
or 

tar xvfz ddscat 6 . l_f ull . tgz 
which wiir'extract" the files from the gzipped tarfile. If your version of "tar" doesn't support the "z" 
option (e.g., you are running Solaris) then try 

zcat ddscat 6 . l_plain . tgz | tar xvf - 
or 

zcat ddscat 6 . l_f ull . tgz | tar xvf - 
If neither of the above work on your system, try the two-stage procedure 
gunzip ddscat 6 . l_plain . tgz 
tar xvf ddscat 6 . l_plain . tar 
or 

gunzip ddscat 6 . l_f ull . tgz 

tar xvf ddscat 6 . l_f ull . tar 
A disadvantage of this two-stage procedure is that it uses more disk space, since after the second step 
you will have the uncompressed tarfile ddscat 6 . l_f ull .tar - about 3.8 Mbytes - in addition to all 
the files you have extracted from the tarfile - another 4.6 Mbytes. 

Any of the above approaches should create subdirectories src, doc, misc, and IDL. The 
source code will be in subdirectory src, and documentation in subdirectory doc. 

5.5 Installing the source code on a MS Windows system 

If you are running Windows on a PC, you will need the "win zip" program]^ win zip should be able to 
"unzip" the gzipped tarfile ddscat 6 . l_plain . tgz or ddscat 6 . l_f ull . tgz and "extract" the 
various files from it automatically. 



^ winzip can be downloaded from |http : / /www . winzip . com| 
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6 Compiling and Linking 
6.1 The "plain" distribution 

In the discussion below, it is assumed that the source code for the "plain" distribution of DDSCAT 
6.1 has been installed in a directory DDA/src. To compile the code on a Linux system with the g77 
compiler, position yourself in the directory DDA/src, and type 
make ddscat 

If you are on another type of system, or wish to use a different compiler (or compiler options) you will 
need to edit the file Makefile to change the choice of compiler (variable FC) compilation options 
(variables FFLAGS and FFLAGSslamcl, and loader options (variable LDFLAGS). 

6.1.1 Compiling the "plain" version on non-Unix systems 

DDSCAT 6.1is written in standard Fortran-77 plus the DO ... ENDDO extension which appears to 
be supported by all current Fortran-77 compatible compilers. It is possible to run DDSCAT on non-Unix 
systems. If the Unix "make" utility is not available, here in brief is what needs to be accomplished: 

All of the necessary Fortran code to compile and link DDSCAT 6.1was included in and should have 
been extracted from the gzipped tarfile ddscat 6 . l_plain . tgz. A subset of the files needs to be 
compiled and linked to produce the ddscat executable. 

The main program is in DDSCAT . f ; it calls a number of subroutines. 

A number of different versions of SUBROUTINE TIME IT are included in ddscat6 . l_plain . tgz. 
The easiest course is to use the version in timeit_null . f , which begins 

SUBROUTINE TIMEIT (CMSGTM, DTIME) 

C 

C timeit_null 
C 

C This version of timeit is a dummy which does not provide any 
C timing information. 

This version does not use any system calls, and therefore the code should compile and link without 
problems; however, when you run the code, it will not report any timing information reporting how 
much time was spent on different parts of the calculation. 

If you wish to obtain timing information, you will need to find out what system calls are supported 
by your operating system. You can look at the other versions of SUBROUTINE TIMEIT provided in 
ddscat 6.1 -plain. tgz to see how this has been done under VMS and various versions of Unix. 

Once you have selected the appropriate version of TIMEIT , you can simply compile and link just as 
you would with any other Fortran code with a number of modules. You should end up with an executable 
with a (system-dependent) name Hke DDSCAT .EXE. 

6.1.2 Optimization 

The performance of DDSCAT 6.1 will benefit from optimization during compilation and the user should 
enable the appropriate compiler flags. 

There is one exception: the LAPACK routine SLAMCl, which is called to determine the number of 
bits of precision provided by the computer architectures, should NOT be optimized, because under some 
optimizers the resulting code will end up in an endless loop. Therefore the Makefile has a separate com- 
pilation flag FFLAGSslamcl which should «of specify optimization (e.g., just g77 -c, or pgf77 
-c).[l 



''Note that the LAPACK routines are only used for some target geometries, and even when used the computation time is insignif- 
icant. 
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6.2 The complete distribution 

In the discussion below, it will be assumed that the complete source code for DDSCAT 6.1 from 
ddscat6 . l_full .tgz has been installed in a directory DDA/ src. To compile the code on a Unix 
system, position yourself in the directory DDA/ src. The as-supplied Make file has compiler options 
set as appropriate for 

• use of the g7 7 compiler under the Linux operating system. 

• use of a Linux-compatible and Solaris-compatible timing routine. 

• FFTW option not enabled 

• MPI not enabled 

• netCDF not enabled 

If you type 

make ddscat 

you should be able to compile and link to create an executable ddscat located in this directory. If 
this does not work with the as-supplied Makefile, the problem is likely due to unavailability of the g7 7 
fortran compiler, or the need to select a different version of the T IME I T routine (see below). If you have 
a different Fortran compiler, you will probably need to edit Makefile to provide the desired compiler 
options. If you are using an operating system other than Linux, you may need to change the Makefile 
choice of timeit. Makefile contains samples of compiler options for selected operating systems, 
including Linux, HP AUX, IBM AIX, and SGI IRIX operating systems. If one of these corresponds to 
your system, try "uncommenting" (i.e., removing the # at the beginning of a line) the appropriate lines 
in the Makefile. 

So far as we know, there are only two operating-system-dependent aspects of DDSCAT 6.1: (1) the 
device number to use for "standard output", and (2) the TIMEIT routine. There are, in addition, three 
installation-dependent aspects to the code: 

• the procedure for linking to the FFTW library (see i]6.2.5b 

• the procedure for linking to the netCDF library (see m0.3\ for discussion of the possibility of 
writing binary files using the machine-independent netCDF standard). 

• the procedure for linking to the MP I library (see ^ 124. 21 ). 

6.2.1 Device Numbers IDVOUT and IDVERR 

The variables IDVOUT and IDVERR specify device numbers for "running output" and "error messages", 
respectively. Normally these would both be set to the device number for "standard output" (e.g., writing 
to the screen if running interactively). Variables IDVERR are set by DATA statements in the "main" 
program DDSCAT . f and in the output routine WRIMSG (file wrimsg . f). The executable statement 
IDVOUT=0 initiahzes IDVOUT to 0. In the as-distributed version of DDSCAT . f , the statement 

OPEN (UNIT=IDVOUT, FILE=CFLLOG) 
causes the output to UNIT=IDVOUT to be buffered and written to the file ddscat . log_nnn, where 
nnn=000 for the first cpu, 001 for the second cpu, etc. If it is desirable to have this output unbuffered 
for debugging purposes, (so that the output will contain up-to-the-moment information) simply comment 
out this OPEN statement. 

6.2.2 Subroutine T IME I T 

The only other operating system-dependent part of DDSCAT 6.1 is the single subroutine TIMEIT. 
Several versions of TIMEIT are provided: 

• timeit_sun . f uses the SunOS system call etime 

• timeit_convex . f uses the Convex OS system call etime 

• timeit_cray . f uses the system call second 
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• timeit_hp . f uses the HP-AUX system calls sysconf and times 

• timeit_ibm6000 . f uses the AIX system call mclock 

• timeit_osf . f uses the DEC OSF system call etime 

• timeit_sgi . f uses the IRIX system call etime 

• timeit_vms . f uses the VMS system calls LIB$ INIT_TIMER and LIB$SHOW_TIMER 

• timeit_titan . f uses the system call cputim 

• timeit_null.fisa dummy routine which provides no timing information, but can be used 
under any operating system. 

You must compile and Unk one of the t ime i t jcxx . f routines, possibly after modifying it to work on 
your system; timeit_null . f is the easiest choice, but it will return no timing information^ 

6.2.3 Optimization 

The performance of DDSCAT 6.1 will benefit from optimization during compilation and the user should 
enable the appropriate compiler flags (e.g., g77 -c -0, or pgf77 -c -fast). 

There is one exception: the LAPACK routine SLAMCl, which is called to determine the number of 
bits of precision provided by the computer architectures, should NOT be optimized, because under some 
optimizers the resulting code will end up in an endless loop. Therefore the Makefile has a separate com- 
pilation flag FFLAGSslamcl which should «of specify optimization (e.g., just g77 -c, or pgf77 
-c).0 

6.2.4 Leaving FFTW Capability Disabled 

FFT1A0 - "Fastest Fourier Transform in the West" - is a publicly available FFT implementation that 
performs very well (see i |T3T l. DDSCAT 6.1 supports use of FFTW, but it is recommended that first-time 
users of DDSCAT first get the code up and running without FFTW support, using the built-in GPFA FFT 
routine written by Clive Temperton, which performs almost as well as FFTW. This is also the option 
you will have to follow if FFTW is not installed on your system (though you can install it yourself, if 
necessary - see i j6.2.5l below). 

The Makefile supplied with the DDSCAT 6.1 distribution is initially set up to leave support for FFTW 
disabled by using a "dummy" subroutine cxf f tw_f ake . f . When you type make ddscat you will 
create a version of ddscat which will not recognize option FFTW2 1 - you mwsf specify option GFPAFT 
in ddscat . par when running the code. 



6.2.5 Enabling FFTW Capability 

The FFTW ("Fastest Fourier Transform in the West") package of Matteo Frigo and Steven G. Johnson 
appears to be one of the fastest public-domain FFT packages available (see i pJl for a performance com- 
parison), and DDSCAT 6.1 allows this package to be used. In order to support the FFTW21 option, 
however, it is necessary to first install the FFTW package on your system (FFTW 2.1.x - we have not 
yet implemented support for FFTW 3.0, which has a different interface). 
FFTW offers two advantages: 

• It appears, for some cases, to be somewhat faster than the GPFA FFT algorithm which is already 
included in the DDSCAT package (see SlT3]and Fig.|6ll. 

Non-Unix sites: The VMS-compatible version of TIMEIT is included (timeit.vms . f). For non-VMS sites, you will 
need to replace this version of TIMEIT with one of the other versions, which are included in ddscat6 . l.plain . tgz and 
ddscat 5 . l_f ull .tgz. If you are having trouble getting this to work, choose the "timeit_null.f ' version of TIMEIT: this will 
return no timing information, but is platform-independent. 

'Note that the LAPACK routines are only used for some target geometries, and even when used the computation time is insignif- 
icant. 

ihttp:/7 www . f f tw . org! 
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• It does not restrict the dimensions of the "computational volume" to be factorizable as 2*3^5*^, 
and therefore for some targets will use a smaller computational volume (and therefore require less 
memory). 

If the FFTW package is not yet installed on your system, then 

1. Download the FFTW 2.1.x package from | http : / / www . f f tw . org| (as of this writing the 
latest version 2.1.x is f f tw-2 .1.5. tar . gz) into some convenient directory. 

2. tar xvfz f f tw-2 . 1 . 5 . tar . gz 

3. cd fftw-2.1.5 

4. ./configure — enable-float — enable-type-pref ix — prefix=PATH 
where PATH is the path to the directory in which you would like the fftw library installed (you 
must have write permission, of course). If you choose to omit 

— pref ix=PATH 

then fftw will be installed in the default directory /usr/local/, but unless you are root you 
may not have write permission to this directory. 

5. make 

6. make install 

This should install the file libsfftw.a in the directory /usr/local/lib (or PATH/lib if you 
used option — pref ix=PATH when running configure). 

Once FFTW is installed, you will need to edit Makefile (go to section 4, "FFTW support": 

• change the lines 

CXFFTW = cxfftw_fake 

LIBFFTW 

to 

#CXFFTW = cxfftw_fake 

#LIBFFTW 

• Go down a few lines in the Makefile, and change the lines 

#CXFFTW = cxfftw 

#LIBFFTW = -L/usr/local/lib -Isfftw 

to 

CXFFTW = cxfftw 

LIBFFTW = -L/usr/local/lib -Isfftw 

LIBFFTW gives the path to libsf ftw . a, which is installation-dependent. Consult your local 
systems administrator. 

Once the above is done, typing make ddscat should create an executable ddscat which will 
recognize option FFTW21 in the file ddscat .par. Note that you will be linking fortran to C, and 
different compilers may behave idiosyncratically. The Makefile has examples which work for g 7 7 and 
pgf 7 7. If you have trouble, consult with someone familiar with linking fortran and C modules on your 
system. 

6.2.6 Leaving netCDF Capability Disabled 

netCDE0is a standard for data transport used at many sites (see SjT03]f or more information on net CDF). 
The Makefile supplied with the distribution of DDSCAT 6.1 is setup to link to a "dummy" subroutine 
writenet_f ake . f instead of subroutine writenet . f , in order to minimize problems during initial 
compilation and linking. The "dummy" routine has no functionality, other than bailing out with a fatal 



http://www.unidata.ucasr.edu/packages/netcdf/ 
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error message if the user makes the mistake of trying to specify one of the netCDF options (ALLCDF or 
ORICDF). First-time users of DDSCAT 6.1 should not try to use the netCDF option - simply skip this 
section, and specify option NOTCDF in ddscat . par. After successfully using DDSCAT 6.1 you can 
return to %.2.7l to try to enable the netCDF capability. 

6.2.7 Enabling netCDF Capability 

Subroutine WRITENET (file writenet . f) provides the capability to output binary data in the netCDF 
standard format (see ijl0.3t . In order to use this routine (instead of writenet_f ake . f), it is necessary 
to take the following steps H 

1 . Have netCDF already installed on your system (check with your system administrator). 

2. Find out where netcdf . inc is located and edit the include statement in writenet . f to 
specify the correct path to netcdf. inc. 

3. Find out where the libnetcdf . a library is located on your system. 

4. Go to section 6 of the Makefile ("netCDF library support") and 

• edit the lines 

WRITENET = writenet_f ake 

LIBNETCDF 

LINKNETCDF 

to 

#WRITENET = writenet_f ake 

#LIBNETCDF 

#LINKNETCDF 

• Go down a few lines in the Makefile and change 

#WRITENET = writenet 

#LIBNETCDF = -L/usr/local/lib 

#LINKNETCDF = -Inetcdf -Insl 
to 

WRITENET = writenet 

LIBNETCDF = -L/usr/local/lib 

LINKNETCDF = -Inetcdf -Insl 

The second and third lines above work on some systems, but will be installation-dependent; if 
you get an error message during the link step, consult your systems administrator or someone 
who has used netCDF on your system. 

6.2.8 Compiling and Linking... 

After suitably editing the Makefile (while still positioned in DDA/src) simply typ^ 
make ddscat 

which should create an executable file DDA/src/ddscat . If the default (as-provided) Makefile is 
used, the g7 7 fortran compiler will be used to create an executable which: 

• will not have MP I support; 

• will not have FFTW support; 

• will not have netCDF capability; 

• will contain timing instructions compatible with Linux, Solaris, as well as several other flavors of 
Unix. 



* Non-UNIX sites: You need to compile the real version of SUBROUTINE WRITENET in writenet . f and not the with the 
fake version in writenet_f ake . f . You will also need to consult your systems administrator to verify that the netCDF library 
has already been installed on your system, and to find out how to link to the library routines. 

' Non-UNIX sites: see ^6^9] 
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To add FFTW capability, see ij6.2.5l To add netCDF capability, see ^6.2.^\ To add MP I support, see 
^241 To replace the Linux-compatible and Sun-compatible timing routine with another, see ji6.2.2l 

6.2.9 Installation of the "full" version on non-Unix systems 

DDSCAT 6.1is written in standard Fortran-77 plus the DO ... ENDDO extension which appears to 
be supported by all current Fortran-77 compatible compilers. It is possible to run DDSCAT on non-Unix 
systems. If the Unix "make" utility is not available, here in brief is what needs to be accomplished: 

All of the necessary Fortran code to compile and link DDSCAT 6.1was included in and should have 
been extracted from the gzipped tarfile ddscat6 . l_full .tgz. A subset of the files needs to be 
compiled and linked to produce the ddscat executable. 

The main program is in DDSCAT . f ; it calls a number of subroutines. 

Some of the subroutines exist in more than one version. 

• Select the version of SUBROUTINE CXFFTW in file cxf ftw_f ake . f instead of the version in 
cxf f tw . f . The resulting code will not support FFTW. 

If you later wish to add FFTW capability, you will first need to install the FFTW library on your 
system (see ! j6.2.5b . After you have done so, you can switch to using the version of SUBROUTINE 
CXFFTW in cxf fwt . f. 

• Select the version of SUBROUTINE WRITENET in writenet.f ake . f instead of the version 
in writenet . f . The resulting code wiU not provide netCDF capability. 

If netCDF capability is required, you will need to install netCDF libraries on your system - see 
%.2Ji . After doing so, you can switch to using writenet . f instead of writenet_f ake . f . 

• Use the code in the file mpi_f ake . f instead of that in the file mpi_f ake . f . The resulting code 
will not support MPl. 

If MPl support is required, you will need to install MPl on your system - see ! j24l After doing so, 
you can switch to using mpi . f . 

• There are a number of different versions of SUBROUTINE TIME IT included in ddscat 6 . l_full .tgz. 
Use the version in timeit_null . f , which begins 

SUBROUTINE TIMEIT (CMSGTM, DTIME) 

C 

C timeit_null 
C 

C This version of timeit is a dummy which does not provide any 
C timing information. 

This version does not use any system calls, and therefore the code should compile and link without 
problems; however, when you run the code, it will not report any timing information reporting how 
much time was spent on different parts of the calculation. 

If you wish to obtain timing information, you will need to find out what system calls are supported 
by your operating system. You can look at the other versions of SUBROUT INE T IMF I T provided 
in ddscat 6 . l_full .tgz to see how this has been done under VMS and various version of 
Unix. 

Once you have selected the appropriate versions of CXFFTW, WRITENET, TIMEIT, and either mpi_f ake . f 
or mpi . f , you can simply compile and link just as you would with any other Fortran code with a number 
of modules. You should end up with an executable with a (system-dependent) name like DDSCAT . EXE. 

In addition to program DDSCAT, we provide one other Fortran 77 program which may be useful. 
Program CALLTARGET (in CALLTARGET . f ) can be used to call the target generation routines. This is 
especially helpful when experimenting with new target geometries. 
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7 Moving the Executable 

Now reposition yourself into the directory DDA (e.g., type cd . .), and copy the executable from 
src/ddscat to the DDA directory by typing 

cp src/ddscat ddscat 
This should copy the file DDA/src/ddscat to DDA/ddscat (you could, of course, simply create a 
symbolic link instead). Similarly, copy the sample parameter file dds cat .par and the file die 1 . t ab 
to the DDA directory by typing 

cp doc/ddscat . par ddscat. par 

cp doc/diel.tab diel.tab 

8 The Parameter File ddscat. par 

The directory DDA should now contain a sample file ddscat . par which provides parameters to the 
program ddscat. As provided (see AppendixlAli. the file ddscat . par is set up to calculate scattering 
by a 32x24 X 16 rectangular array of 12288 dipoles, with an effective radius a^ff ~ 2fj,m, at a wavelength 
of 6.2832/im (for a "size parameter" 27racff/A = 2). 

The dielectric function of the target material is provided in the file diel . tab, which is a sample 
file in which the refractive index is set to m = 1.33 + O.Olz at all wavelengths; the name of this file is 
provided to ddscat by the parameter file ddscat . par. 

The sample parameter file as supplied calls for the GPFA FFT routine (GPFAFT) of Temperton 
(1992) to be employed and the PBCGST iterative method to be used for solving the system of linear 
equations. (See section fJT3]and SjT4]for discussion of choice of FFT algorithm and choice of equation- 
solving algorithm.) 

The sample parameter file specifies (via option LATTDR) that the "Lattice Dispersion Relation" 
(Draine & Goodman 1993) be employed to determine the dipole polarizabilities. See ifTTI for discussion 
of other options. 

The sample parameter files specifies options NOTBIN and NOTCDF so that no binary data files and 
no NetCDF files will be created by ddscat. The sample ddscat .par file specifies that the 
calculations be done for a single wavelength (6.2832/xm) and a single effective radius (2.00/im). Note 
that in DDSCAT 6.1 the "effective radius" Ooft is the radius of a sphere of equal volume - i.e., a sphere 
of volume NcP , where d is the lattice spacing and N is the number of occupied (i.e., non-vacuum) 
lattice sites in the target. Thus the effective radius Ocff = {3N/4ny/^d . 

The incident radiation is always assumed to propagate along the x axis in the "Lab Frame". The sam- 
ple ddscat . par file specifies incident polarization state Bqi to be along the y axis (and consequently 
polarization state §02 will automatically be taken to be along the z axis). I0RTH=2 in ddscat .par 
calls for calculations to be carried out for both incident polarization states (gqi and eo2 - see i|2Tl i. 

The target is assumed to have two vectors ai and 3.2 embedded in it; 3.2 is perpendicular to ai . In 
the case of the 32x24x 16 rectangular array of the sample calculation, the vector ai is along the "long" 
axis of the target, and the vector 3.2 is along the "intermediate" axis. The target orientation in the Lab 
Frame is set by three angles: /?, Q, and $, defined and discussed below in SjlT] Briefly, the polar angles 
G and $ specify the direction of ai in the Lab Frame. The target is assumed to be rotated around ai 
by an angle (3. The sample ddscat .par file specifies (3 = and $ = (see lines in ddscat .par 
specifying variables BETA and PHI), and calls for three values of the angle Q (see line in ddscat . par 
specifying variable THETA). DDSCAT 6.1 chooses Q values uniformly spaced in cos 6; thus, asking 
for three values of 9 between and 90° yields 9 = 0, 60°, and 90°. 

Appendix lAlprovides a detailed description of the file ddscat . parFl 



Note that the structure of ddscat . par depends on the version of DDSCAT being used. Make sure you update old parameter 
files before using them with DDSCAT 6.1 ! 
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9 Running DDSCAT 6.1 Using the Sample dds cat . par File 

To execute the program on a UNIX system (running either sh or csh), simply type 

ddscat >& ddscat.out & 
which will redirect the "standard output" to the file ddscat . out, and run the calculation in the back- 
ground. The sample calculation [32x24x16=12288 dipole target, 3 orientations, two incident polar- 
izations, with scattering (Mueller matrix elements Sij) calculated for 38 distinct scattering directions], 
requires 27.9 cpu seconds to complete on a 2.0 GHz Xeon workstation. 

10 Output Files 
10.1 ASCII files 

If you run DDSCAT using the command 
ddscat >& ddscat.out & 
you will have various types of ASCII files when the computation is complete: 

• a file ddscat . out; 

• a file ddscat . log_000; 

• afilemtable; 

• a file qtable; 

• a file qtable2; 

• files vxxxryyori . avg (one, wOOOrOOori . avg, for the sample calculation); 

• if ddscat .par specified IWRKSC=1, there will also be files vxxxryykzzz . sea (3 for the sample 
calculation: wOOOrOOkOOO . sea, wOOOrOOkOOl . sea, w000r00k002 . sea). 

The file ddscat . out will contain minimal information (it may in fact be empty). 

The file ddscat . log_000 will contain any error messages generated as well as a running report 
on the progress of the calculation, including creation of the target dipole array. During the iterative 
calculations, Qcxt, Qabs, and Qpha are printed after each iteration; you will be able to judge the degree 
to which convergence has been achieved. Unless TIME IT has been disabled, there will also be timing 
information. If the MP I option is used to run the code on multiple cpus, there will be one file of the form 
ddscat . log_nnn for each of the cpus, with nnn=000 , 001,002, . . .. 

The file mtable contains a summary of the dielectric constant used in the calculations. 

The file qtable contains a summary of the orientationally-averaged values of Qcxt, Qabs, Qsca, 
g{l) = {cos{9s)), (cos^(0s)), Qbk, and N^ca- Here Qext, Qabs, and Qsca are the extinction, absorption, 
and scattering cross sections divided by Tra^jj. Qbk is the differential cross section for backscattering 
(area per sr) divided by Tra^g. Agca is the number of scattering directions used for averaging over 
scattering directions (to obtain (cos 6), etc.) (see ^2% . 

The file qtable2 contains a summary of the orientationally-averaged values of Qpha, Qpoi, and 
Qcpoi- Here Qpha is the "phase shift" cross section divided by Tra^g (see definition in Draine 1988). 
Qpoi is the "polarization efficiency factor", equal to the difference between Qcxt for the two orthogonal 
polarization states. We define a "circular polarization efficiency factor" Qcpoi = QpoiQpha, since an 
optically-thin medium with a small twist in the alignment direction will produce circular polarization in 
initially unpolarized light in proportion to Qcpoi- 

For each wavelength and size, DDSCAT 6.1 produces a file with a name of the form 
wxxxryyori . avg, where index xxx (=000, 001, 002....) designates the wavelength and index yy (=00, 
01, 02...) designates the "radius"; this file contains Q values and scattering information averaged over 
however many target orientations have been specified (see fJTT] The file wOOOrOOori.avg produced 
by the sample calculation is provided below in AppendixlB] 

In addition, if ddscat .par has specified IWRKSC=1 (as for the sample calculation), DDSCAT 
6.1 will generate files with names of the form wxxxryykzzz . avg, where xxx and yy are as before, and 
index zzz =(000,001,002...) designates the target orientation; these files contain Q values and scattering 
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information for each of the target orientations. The structure of each of these files is very similar to 
that of the vixxxryyori . avg files. Because these files may not be of particular interest, and take up 
disk space, you may choose to set IWRKSC=0 in future work. However, it is suggested that you run the 
sample calculation with IWRKSC=1. 

The sample ddscat . par file specifies IWRKSC=1 and calls for use of 1 wavelength, 1 target size, 
and averaging over 3 target orientations. Running DDSCAT 6.1 with the sample ddscat . par file will 
therefore generate files wOOOrOOkOOO . sea, wOOOrOOkOOl . sea, and w000r00k002 . sea . To 
understand the information contained in one of these files, please consult Appendix |C] which contains 
an example of the file wOOOrOOkOOO.sca produced in the sample calculation. 

10.2 Binary Option 

It is possible to output an "unformatted" or "binary" file (dd.bin) with fairly complete information, 
including header and data sections. This is accomplished by specifying either ALLBIN or ORIBIN in 
ddscat . par . 

Subroutine writebin . f provides an example of how this can be done. The "header" section con- 
tains dimensioning and other variables which do not change with wavelength, particle geometry, and tar- 
get orientation. The header section contains data defining the particle shape, wavelengths, particle sizes, 
and target orientations. If ALLBIN has been specified, the "data" section contains, for each orientation, 
Mueller matrix results for each scattering direction. The data output is limited to actual dimensions of 
arrays; e.g. nscat ,4,4 elements of Muller matrix are written rather than mxscat ,4,4. This is an 
important consideration when writing postprocessing codes. 

A skeletal example of a postprocessing code was written by us (PJF) and is provided in subdirectory 
DDA/IDL. If you do plan to use the Interactive Data Language (IDL) for postprocessing, you may 
consider using the netCDF binary file option which offers substantial advantages over the FORTRAN 
unformatted write. More information about IDL is available at http : / /www . rsinc . com/idl. 
Unfortunately IDL requires a license and hence is not distributed with DDSCAT. 

10.3 Machine-Independent Binary File Option: netCDF 

The "unformatted" binary file is compact, fairly complete, and (in comparison to ASCII output files) is 
efficiently written from FORTRAN. However, binary files are not compatible between different computer 
architectures if written by "unformatted" WRITE from FORTRAN. Thus, you have to process the data 
on the same computer architecture that was used for the DDSCAT calculations. We provide an option of 
using netCDF with DDSCAT - if the user specifies either the ALLCDF or ORICDF options, DDSCAT 
will write an output file with output written to a file (dd . cdf ) in netCDF format (in addition to the usual 
ASCII output files). 

The netCDF librarjji^ defines a machine-independent format for representing scientific data. To- 
gether, the interface, library, and format support the creation, access, and sharing of scientific data. For 
more information, go to the net cdf websitJUl 

Several major graphics packages (for example IDL) have adopted netCDF as a standard for data 
transport. For this reason, and because of strong and free support of netCDF over the network by 
UNIDATA, we have implemented a netCDF interface in DDSCAT. This may become the method of 
choice for binary file storage of output from DDSCAT. 

After the initial "learning curve" netCDF offers many advantages: 

• It is easy to examine the contents of the file (using tools such as ncdump). 

• Binary files are portable - they can be created on a supercomputer and processed on a workstation. 

• Major graphics packages now offer netCDF interfaces. 

• Data input and output are an order of magnitude faster than for ASCII files. 

• Output data files are compact. 



" fhttp r/ / www . unidata . ucar . edu /packages /net cdf /| 
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X = 2Tra/X 



Figure 5: Scattering and absorption for a m = 1.7 + O.li sphere, calculated using two prescriptions for the polarizability: 
LATTDR is the lattice dispersion relation result of Draine & Goodman (1993). GKDLDR is the lattice dispersion relation result of 
Gutkowicz-Krusin & Draine (2004). We see that there is little practical difference between them for this case. For other examples 
see Gutkowicz-Krusin & Draine (2004). 



The disadvantages include: 

• Need to have netCDF installed on your system. 

• Lack of support of complex numbers. 

• Nontrivial learning curve for using netCDF. 

The public-domain netCDF software is available for many operating systems from 
|htt p : /"/www . unidata . ucar . edu/packages/netcdf . The steps necessary for enabling the 
netCDF capability in DDSCAT 6.1 are listed above in %.2J\ Once these have been successfully ac- 
compUshed, and you are ready to run DDSCAT to produce netCDF output, you must choose either 
the ALLCDF or ORICDF option in ddscat .par; ALLCDF will result in a file which will include the 
MuUer matrix for every wavelength, particle size, and orientation, whereas ORICDF will result in a file 
limited to just the orientational averages for each wavelength and target size. 



11 Dipole Polarizabilities 

Option LATTDR specifies that the "Lattice Dispersion Relation" of Draine & Goodman (1993) be em- 
ployed to determine the dipole polarizabilities. This polarizability works well. 

Option GKDLDR specifies that the polarizability be prescribed by the "Lattice Dispersion Relation", 
but with the polarizability found by Gutkowicz-Krusin & Draine (2004), who corrected an error in the 
analysis of Draine & Goodman (1993). The GKDLDR polarizabiHty differs somewhat from the LATTDR 
polarizability, but the differences in calculated scattering cross sections are relatively small, as can be 
seen from Figure |5] 



DIELECTRIC FUNCTIONS 



12 Dielectric Functions 

12.1 Table Specifying Dielectric Function or Refractive Index 

In order to assign the appropriate dipole polarizabilities, DDSCAT 6.1 must be given the dielectric 
constant of the material (or materials) of which the target of interest is composed. Unless the user wishes 
to use the dielectric function of either solid or liquid H2O (see below), this information is supplied to 
DDSCAT through a table (or tables), read by subroutine DIELEC in file dielec . f , and providing 
either the complex refractive index m = n + ik 01 complex dielectric function e = ei+ ie2 as a function 
of wavelength A. Since m = e^/^, or e = m^, the user must supply either m or e, but not both. 

The table formatting is intended to be quite flexible. The first line of the table consists of text, up 
to 80 characters of which will be read and included in the output to identify the choice of dielectric 
function. (For the sample problem, it consists of simply the statement m = 1.33 + O.Oli.) The 
second line consists of 5 integers; either the second and third or the fourth and fifth should be zero. 

• The first integer specifies which column the wavelength is stored in. 

• The second integer specifies which column Re(m) is stored in. 

• The third integer specifies which column Im(TO) is stored in. 

• The fourth integer specifies which column Re(e) is stored in. 

• The fifth integer specifies which column Im(e) is stored in. 

If the second and third integers are zeros, then DIELEC will read Re(e) and Im(e) from the file; if the 
fourth and fifth integers are zeros, then Re(m) and lm{77i) will be read from the file. 

The third line of the file is used for column headers, and the data begins in line 4. There must be at 
least 3 lines of data: even if e is required at only one wavelength, please supply two additional "dummy" 
wavelength entries in the table so that the interpolation apparatus will not be confused. 

12.2 Multiple Dielectric Functions 

For targets composed of more than one material, or composed of anisotropic materials, it is necessary to 
specify more than one dielectric function. This is accomplished by providing multiple files; one file per 
dielectric function. The files are identified in ddscat . par. Here is an exa mple f or a target consisting 
of two concentric ellipsoids (see description of target option CONELL in iil9.3l ): the inner ellipsoid 
consists of material with dielectric function (or refractive index) tabulated in the file diel_l .tab, 
while the region between the inner ellipsoid and outer ellipsoid has dielectric function (or refractive 
index) tabulated in the file diel_2 . tab. Note that ddscat . par specifies NCOMP (the number of 
dielectric materials to be 2. 

' =================== Parameter file ===================' 

PRELIMINARIES ****' 

'NOTORQ' = CMT0RQ*6 (DOTORQ, NOTORQ) — either do or skip torque calculations 

'PBCGST' = CMDS0L*6 (PBCGST, PETRKP) — select solution method 

'GPFAFT' = CMETHD*6 (GPFAFT, FFTWJ, CONVEX) 

'LATTDR' = CALPHA*6 (LATTDR, SCLDR) 

'NOTBIN' = CBINFLAG (ALLBIN, ORIBIN, NOTBIN) 

'NOTCDF' = CNETFLAG (ALLCDF, ORICDF, NOTCDF) 

'CONELL' = CSHAPE*6 ( FRMF IL , ELLIPS , C YLNDR, RCTNGL , HEXGON , TETRAH , UNIC YL, . . . ) 
32 24 16 16 12 8 = shape parameters SHPARl, SHPAR2, SHPAR3, . . . 
2 = NCOMP = number of dielectric materials 

'TABLES' = CDIEL*6 ( TABLES , H20ICE, H20LIQ ; if TABLES, then filenames follow...) 
'diel_l.tab' = name of file containing dielectric function 
' diel_2 .tab' 

'**** CONJUGATE GRADIENT DEFINITIONS ****' 
= INIT (TO BEGIN WITH | X0> = 0) 

l.OOe-5 = TOL = MAX ALLOWED (NORM OF | G>=AC I E>-ACA | X> ) / (NORM OF AC|E>) 
' **** Angular resolution for calculation of <cos>, etc. +***' 
0.5 = ETASCA (number of angles is proportional to [ ( 3+x) /ETASCA] " 2 ) 

' **** Wavelengths (micron) ****' 

6.283185 6.283185 1 ' INV = wavelengths (first, last, how many, how=LIN, INV, LOG) 
' **** Effective Radii (micron) **** ' 

2 2 1 'LIN' = eff. radii (first, last, how many, how=LIN, INV, LOG) 
**** Define Incident Polarizations *+++' 
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(0,0) (l.,0.) (0.,0.) = Polarization state eOl (k along x axis) 

2 = lORTH (=1 to do only pol . state eOl; =2 to also do orth. pol . state) 

I = IWRKSC (-0 to suppress, =1 to write ".sea" file for each target orient. 
'**** Prescribe Target Rotations ****' 

0. 0. 1 = BETAMI, BETAMX, NBETA (beta=rotation around al) 

0. 90. 3 = THETMI, THETMX, NTHETA (theta=angle between al and k) 

0. 0. 1 = PHIMIN, PHIMAX, NPHI (phi=rotation angle of al around k) 

Specify first IWAV, IRAD, lORI (normally 0) 
= first IWAV, first IRAD, first lORI (0 to begin fresh) 
' Select Elements of S_ij Matrix to Print 

6 = NSMELTS = number of elements of S_ij to print (not more than 9) 

II 12 21 22 31 41 = indices ij of elements to print 

Specify Scattered Directions ****' 
0. 0. 180. 10 = phi, thetan_min, thetan_max, dtheta (in degrees) for plane A 
90. 0. 180. 10 = phi, ... for plane B 



12.3 Water Ice or Liquid Water 

In the event that the user wishes to compute scattering by targets with the refractive index of either solid 
or liquid H2O, we have included two "built-in" dielectric functions. If H20ICE is specified on line 10 of 
dds cat . par, DDSCAT will use the dielectric function of H2O ice at T=250K as compiled by Warren 
(1984). If H20LIQ is specified on line 10 of ddscat .par, DDSCAT will use the dielectric function 
for liquid H2O at T=280K using subroutine REFWAT written by Eric A. Smith. For more information 
see |http : / / atol . ucsd . edu/^pf latau/ re frtab/ Index . htni| . 



13 Choice of FFT Algorithm 

One major change in going from DDSCAT.4b to 4c was modification of the code to permit use of the 
GPFA FFT algorithm developed by Dr Clive Temperton (1992). In going from DDSCAT 5al0 to 6.0 
we introduced a new FFT option: the FFTW package of Frigo and Johnson. DDSCAT 6.1continues to 
offer Temperton's GPFA code as an alternative (and quite fast) FFT implementation. 

Table [T] is the res . dat file created when run on a 2.0 GHz Intel Xeon system, using the pgf77 
compiler: 

It is clear that both the GPFA code and the FFTW code are much faster than the Brenner FFT (by 
a factor of 7 for the lOOx lOOx 100 case). The FFTW code and GPFA code are quite comparable in 
performance - for some cases the GPFA code is faster, for other cases the FFTW code is faster. For 
target dimensions which are factorizable as 2'3''5'^ (for integer i, j, k), the GPFA and FFTW codes 
have the same memory requirements. For targets with extents N^, Ny, which are not factorizable 
as 2*3^5'^, the GPFA code needs to "extend" the computational volume to have values of Nx, Ny, and 
Nz which are factorizable by 2, 3, and 5. For these cases, GPFA requires somewhat more memory than 
FFTW. However, the fractional difference in required memory is not large, since integers factorizable as 
2*3^ 5'^' occur fairly frequently[3 

Based on Figure 5, we see that while for some cases FFTW 2 . 1 . 5 is faster than the GPFA algo- 
rithm, the difference appears to be fairly marginal. 

The GPFA code contains a parameter LVR which is set in data statements in the routines gpf a2 f , 
gpf a3f, and gpf a5f . LVR is supposed to be optimized to correspond to the "length of a vector reg- 
ister" on vector machines. As delivered, this parameter is set to 64, which is supposed to be appropriate 
for Grays other than the C90 (for the C90, 128 is supposed to be preferable)13 The value of LVR is not 
critical for scalar machines, as long as it is fairly large. We found little difference between LVR=64 and 
128 on a Sparc 10/51, on an Ultrasparc 170, and on an Intel Xeon cpu. You may wish to experiment 

'^2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 40, 45, 48, 50, 54, 60, 64, 72, 75, 80, 81, 90, 96, 100, 108, 120, 
125, 128, 135, 144, 150, 160, 162, 180, 192, 200 are the integers < 200 which are of the form 2'3^ 5'=. 

In place of "preferable" users are encouraged to read "necessary" - there is some basis for fearing that results computed on a 
C90 with LVR other than 128 run the risk of being incorrect! Please use LVR=1 2 8 if running on a C90! 
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Table 1: Timings for 3 different 3-d FFT implementations 

CPU time (sec) for different 3-D FFT methods 
Machine=2 . GHz Xeon, 256k:B L2 cache, pgf77 compiler 
parameter LVR = 64 

LVR="length of vector registers" in gpfa2f , gpf a3f , gpf a5f 



Brenner Temperton Frigo S Johnson 



NX 


NY 


NZ 


(FOURX) 


(GPFA) 


(FFTW) 


31 


31 


31 


0, 


.040 


, 


. 082 


0, 


.030 


32 


32 


32 


0, 


. 012 


, 


. 007 


0, 


.004 


34 


34 


34 


0, 


.040 


, 


. 120 


0, 


.023 


36 


36 


36 


0, 


. 047 


, 


. 008 


0, 


.008 


38 


38 


38 


0, 


.060 


, 


. 013 


0, 


.037 


40 


40 


40 


0, 


. 055 


, 


. 013 


0, 


.008 


43 


43 


43 


0, 


. 155 


, 


. 018 


0, 


. 072 


45 


45 


45 


0, 


. 125 


, 


. 017 


0, 


. 014 


47 


47 


47 


0, 


.220 


, 


. 023 


0, 


. 132 


48 


48 


48 


0, 


. 120 


, 


. 023 


0, 


.016 


49 


49 


49 


0, 


. 150 


, 


. 024 


0, 


.019 


50 


50 


50 


0, 


. 170 


, 


. 024 


0, 


.020 


52 


52 


52 


0, 


.180 


, 


. 033 


0, 


.028 


54 


54 


54 


0, 


.260 


, 


. 031 


0, 


.023 


57 


57 


57 


0, 


.310 


, 


. 042 


0, 


. 147 


60 


60 


60 


0, 


.320 


, 


. 042 


0, 


. 042 


62 


62 


62 


0, 


.440 


, 


.068 


0, 


.210 


64 


64 


64 


0, 


.240 


, 


.068 


0, 


. 055 


68 


68 


68 


0, 


.470 


, 


.090 


0, 


.190 


72 


72 


72 


0, 


. 600 


, 


.093 


0, 


.060 


74 


74 


74 


0, 


.850 


, 


.093 


0, 


. 347 


75 


75 


75 


0, 


.730 


, 


.093 


0, 


. 077 


78 


78 


78 


0, 


.800 


, 


. 120 


0, 


.090 


80 


80 


80 


0, 


.710 


, 


. 120 


0, 


.100 


81 


81 


81 


1 , 


.070 


, 


. 115 


0, 


.100 


86 


86 


86 


1 , 


. 510 


, 


. 170 


0, 


. 630 


90 


90 


90 


1 , 


.390 


, 


. 170 


0, 


. 125 


93 


93 


93 


1 , 


.740 


, 


.260 


0, 


. 720 


96 


96 


96 


1 , 


. 920 


, 


.260 


0, 


.270 


98 


98 


98 


1 , 


.590 


, 


.230 


0, 


. 170 


100 


100 


100 


1 , 


. 670 


, 


.225 


0, 


.225 


104 


104 


104 


1 , 


.740 


, 


.280 


0, 


.220 


108 


108 


108 


2 , 


.490 


, 


.280 


0, 


.210 


114 


114 


114 


2 , 


. 940 


, 
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1 , 


.100 


120 
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3, 


.100 


, 
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0, 
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, 
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.750 
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3, 
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, 


.450 


0, 
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, 
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3, 
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, 
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1 , 
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, 
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0, 
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. 420 


0, 
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0, 
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140 
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0, 
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0, 
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, 
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0, 
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, 
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0, 
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, 
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0, 
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155 
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3, 
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8, 
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3, 
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0, 
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171 
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.380 
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3, 


.770 
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180 
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192 
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.050 
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3, 
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196 
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Figure 6: Comparison of cpu time required by 3 different FFT implementations. It is seen that the GPFA and FPTW implemen- 
tations have comparable speeds, much faster than Brenner's FFT implementation. 

with different LVR values on your computer architecture. To change LVR, you need to edit gpf a . f and 
change the three data statements where LVR is set. 

The choice of FFT implementation is obtained by specifying one of: 

• FFTW21 to use the FFTW 2.1.x algorithm of Frigo and Johnson - this is recommended, but 
requires that the FFTW 2.1.x library be installed on your system. 

• GPF AFT to use the GPFA algorithm (Temperton 1992) - this is almost as fast as the FFTW 
package, and does not require that the FFTW package be installed on your system; 

• CONVEX to use the native FFT routine on a Convex system (this option is presumably obsolete, 
but is retained as an example of how to use a system-specific native routine). 

14 Choice of Iterative Algorithm 

As discussed elsewhere (e.g., Draine 1988), the problem of electromagnetic scattering of an incident 
wave Einc by an array of N point dipoles can be cast in the form 



where E is a 37V-dimensional (complex) vector of the incident electric field Einc at the N lattice sites, 
P is a 3iV-dimensional (complex) vector of the (unknown) dipole polarizations, and A is a 3A^ x 3A^ 
complex matrix. 

Because 3N is a large number, direct methods for solving this system of equations for the unknown 
vector P are impractical, but iterative methods are useful: we begin with a guess (typically, P = 0) for 
the unknown polarization vector, and then iteratively improve the estimate for P until equation (flTl l is 
solved to some error criterion. The error tolerance may be specified as 



AP = E 



(17) 



|AtAP- AtE| 
|AtE| 



< h 



(18) 
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where is the Hermitian conjugate of A [{A'')ij = (^1^;)*], and h is the error tolerance. We typically 
use h = 10^^ in order to satisfy eq.lfTTb to high accuracy. The error tolerance h can be specified by the 
user through the parameter TOL in the parameter file ddscat . par (see AppendixlAl. 

A major change in going from DDSCAT.4b to 5a (and subsequent versions) was the implementation 
of several different algorithms for iterative solution of the system of complex linear equations. DDSCAT 
6.1 is now structured to permit solution algorithms to be treated in a fairly "modular" fashion, facilitating 
the testing of different algorithms. Many algorithms were compared by Flatau (19970; two of them 
performed well and are made available to the user in DDSCAT 6.1. The choice of algorithm is made by 
specifying one of the options: 

• PBCGST - Preconditioned BiConjugate Gradient with STabilitization method from the Parallel 
Iterative Methods (PIM) package created by R. Dias da Cunha and T. Hopkins. 

• PETRKP - the complex conjugate gradient algorithm of Petravic & Kuo-Petravic (1979), as coded 
in the Complex Conjugate Gradient package (CCGPACK) created by P.J. Flatau. This is the algo- 
rithm discussed by Draine (1988) and used in previous versions of DDSCAT. 

Both methods work well. Our experience suggests that PBCGST is often faster than PETRKP, by perhaps 
a factor of two. We therefore recommend it, but include the PETRKP method as an alternative. 

The Parallel Iterative Methods (PIM) by Rudnei Dias da Cunha (rddSukc . ac . uk) and Tim Hop- 
kins (trh@ukc.ac. uk) is a collection of Fortran 77 routines designed to solve systems of linear equa- 
tions on parallel and scalar computers using a variety of iterative methods (available at 
|http : / /www . mat . uf rgs . br/pim-e . html| l. PIM offers a number of iterative methods, includ- 
ing 

• Conjugate-Gradients, CG (Hestenes 1952), 

• Bi-Conjugate-Gradients, BICG (Fletcher 1976), 

• Conjugate-Gradients squared, CGS (Sonneveld 1989), 

• the stabilised version of Bi-Conjugate-Gradients, BICGSTAB (van der Vorst 1991), 

• the restarted version of BICGSTAB, RBICGSTAB (Sleijpen & Fokkema 1992) 

• the restarted generalized minimal residual, RGMRES (Saad 1986), 

• the restarted generalized conjugate residual, RGCR (Eisenstat 1983), 

• the normal equation solvers, CGNR (Hestenes 1952 and CGNE (Craig 1955), 

• the quasi-minimal residual, QMR (highly parallel version due to Bucker & Sauren 1996), 

• transpose-free quasi-minimal residual, TFQMR (Freund 1992), 

• the Chebyshev acceleration, CHEBYSHEV (Young 1981). 

The source code for these methods is distributed with DDSCAT but only PBCGST and PETRKP can 
be called directly via ddscat. par. It is possible to add other options by changing the code in 
get f ml . f . Flatau (1997) has compared the convergence rates of a number of different methods. 

A helpful introduction to conjugate gradient methods is provided by the report "Conjugate Gradient 
Method Without Agonizing Pain" by Jonathan R. Shewchuk0 



15 Calculation of (cos 9), Radiative Force, and Radiation Torque 

In addition to solving the scattering problem for a dipole array, DDSCAT 6.1 can compute the three- 
dimensional force Frad and torque Fiad exerted on this array by the incident and scattered radiation 
fields. The radiation torque calculation is carried out, after solving the scattering problem, only if 
DOTORQ has been specified in ddscat .par. For each incident polarization mode, the results are 
given in terms of dimensionless efficiency vectors Q^^ and Qp, defined by 

Qpr ^ 2 ' (19) 

A postscript copy of this report - file eg . ps - is distributed with the DDSCAT 6.1 documentation. 

Available as a postscript file: [ftp :/ /REPORTS .ADM . CS . CMU . EDU/usr /anon/1 994 /CMU-CS-94-12 5 . ps| 



MEMORY REQUIREMENTS 



Qr^-?^ , (20) 

where Frad and Fiad are the time-averaged force and torque on the dipole array, k = 27r/A is the 
wavenumber in vacuo, and Uiad = Eq/8tt is the time-averaged energy density for an incident plane 
wave with amplitude Eq cos{ujt + (f). The radiation pressure efficiency vector can be written 

Qpr = Qcxtk - Qscag , (21) 

where k is the direction of propagation of the incident radiation, and the vector g is the mean direction 
of propagation of the scattered radiation: 

1 /■,^ rfC.ca(n,k) . 

where dVL is the element of solid angle in scattering direction n, and dCsca/dil is the differential scat- 
tering cross section. 

Equations for the evaluation of the radiative force and torque are derived by Draine & Weingartner 
(1996). It is important to note that evaluation of Qpr and Qr involves averaging over scattering di- 
rections to evaluate the linear and angular momentum transport by the scattered wave. This evaluation 
requires appropriate choices of the parameter ETASCA- see ^22\ 

In addition, DDSCAT 6.1 calculates (cos 9) [the first component of the vector g in eq. (l22l il and the 
second moment (cos^ 9) . These two moments are useful measures of the anisotropy of the scattering. 
For example, Draine (2003) gives an analytic approximation to the scattering phase function of dust 
mixtures that is parameterized by the two moments (cos 6) and (cos^ 9) . 



16 Memory Requirements 

Since Fortran-77 does not allow for dynamic memory allocation, the executable has compiled into it the 
dimensions for a number of arrays; these array dimensions limit the size of the dipole array which the 
code can handle. The source code supplied to you (which can be used to run the sample calculation 
with 12288 dipoles) is restricted to problems with targets with a maximum extent of 32 lattice spacings 
along the x-, y-, and z-directions (MXNX=32 , MXNY=32 , MXNZ = 32; i.e, the target must fit within an 
32x32x32=32768 cube) and involve at most 9 different dielectric functions (MXC0MP=9). With this 
dimensioning, the executable requires about 1.3 MB of memory to run on an Ultrasparc system; memory 
requirements on other hardware and with other compilers should be similar. To run larger problems, you 
will need to edit the DDSCAT . f to change certain PARAMETER values and then recompile. 

All of the dimensioning takes place in the main program DDSCAT - this should be the only routine 
which it is necessary to recompile. All of the dimensional parameters are set in PARAMETER statements 
appearing before the array declarations. You need simply edit the parameter statements. Remember, of 
course, that the amount of memory allocated by the code when it runs will depend upon these dimen- 
sioning parameters, so do not set them to ridiculously large values! The parameters MXNX, MXNY, MXNZ 
specify the maximum extent of the target in the x-, y-, or z-directions. The parameter MXCOMP speci- 
fies the maximum number of different dielectric functions which the code can handle at one time. The 
comment statements in the code supply all the information you should need to change these parameters. 

The memory requirement for DDSCAT 6.1 (with the netCDF and FFTW options disabled) is ap- 
proximately 

(1690 + 0.764MXNX X MXNY X MXNZ) kbytes (23) 

Thus a 32x32x32 calculation requires 22.5 MBytes, while a 48x48x48 calculation requires 73.0 
MBytes. These values are for the pgf 7 7 compiler on an Intel Xeon system running the Linux oper- 
ating system. 
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Fi gure 7 ! Target orientation in the Lab Frame, x is the direction of propagation of the incident radiation, and y is the direction of 
the "real" component of the first incident polarization mode. In this coordinate system, the orientation of target axis ai is specified 
by angles B and $. With target axis ai fixed, the orientation of target axis a.2 is then determined by angle /3 specifying rotation of 
the target around ai. When /3 = 0, 3.2 lies in the ai,x plane. 



17 Target Orientation 

Recall that we define a "Lab Frame" (LF) in which the incident radiation propagates in the +x direction. 
In ddscat . par one specifies the first polarization state eoi (which obviously must lie in the y, z plane 
in the LF); DDSCAT automatically constructs a second polarization state eo2 = x x Gqj^ orthogonal to 
eoi (here x is the unit vector in the +a; direction of the LF. For purposes of discussion we will always let 
unit vectors x, y, z = x x y be the three coordinate axes of the LF. Users will often find it convenient to 
let polarization vectors Gqi = y- = z (although this is not mandatory - see jjSTT i. 

Recall that definition of a target involves specifying two unit vectors, ai and a2, which are imagined 
to be "frozen" into the target. We require a2 to be orthogonal to ai. Therefore we may define a "Target 
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Frame" (TF) defined by the three unit vectors ai, 3.2, and as = ai x a2 . 

For example, when DDSCAT creates a 32x24x16 rectangular solid, it fixes ai to be along the 
longest dimension of the solid, and 3.2 to be along the next-longest dimension. 

Orientation of the target relative to the incident radiation can in principle be determined two ways: 

1 . specifying the direction of ai and a2 in the LF, or 

2. specifying the directions of x (incidence direction) and y in the TF. 

DDSCAT 6.1 uses method 1.: the angles 8, and (3 are specified in the file ddscat . par. The target 
is oriented such that the polar angles Q and $ specify the direction of ai relative to the incident direction 
X, where the x,y plane has $ = 0. Once the direction of ai is specified, the angle f3 than specifies how 
the target is to rotated around the axis ai to fully specify its orientation. A more extended and precise 
explanation follows: 



17.1 Orientation of the Target in the Lab Frame 

DDSCAT uses three angles, 8, <I>, and f3, to specify the directions of unit vectors ai and a2 in the LF 
(see Fig.|7]i. 

8 is the angle between ai and x. 

When $ = 0, ai will lie in the x, y plane. When $ is nonzero, it will refer to the rotation of ai 
around x: e.g., $ = 90° puts ai in the x, z plane. 

When /? = 0, 3.2 will lie in the x, ai plane, in such a way that when 8 = and $ = 0, 3.2 is in the y 
direction: e.g, 8 = 90°, $ = 0, /? = has si = y and a2 = — x. Nonzero /3 introduces an additional 
rotation of a2 around ai: e.g., 8 = 90°, <1> = 0, = 90° has ai = y and a2 = z. 

Mathematically: 

ai = X cos 8 + y sin 8 cos $ + z sin 8 sin $ (24) 

3.2 = —X sin 8 cos f3 + y [cos 8 cos /3 cos $ — sin (3 sin <1>] 

+z[cos 8 cos/3sin $ + sin/3 cos $] (25) 

3.3 = X sin 8 sin /3 — y [cos 8 sin (3 cos $ + cos /3 sin ^] 

— z[cos8sin/3sin<i> — cos/3cos$] (26) 



or, equivalently: 



X = 3i cos 8 — 3.2 sin 8 cos (3 + sls sin 8 sin (3 (27) 
y = 3i sin 8 cos $ + 3.2 [cos 8 cos /3 cos $ — sin (3 sin <i>] 

— a3[cos8sin/3cos$ + cos/?sin$] (28) 
z = 3.1 sin 8 sin $ + a2 [cos 8 cos (3 sin <I> + sin (3 cos $] 

— a3[cos8sin/3sin<l> — cos/3cos$] (29) 



17.2 Orientation of the Incident Beam in the Target Frame 

Under some circumstances, one may wish to specify the target orientation such that x (the direction of 
propagation of the radiation) and y (usually the first polarization direction) and z (= x x y) refer to 
certain directions in the TF. Given the definitions of the LF and TF above, this is simply an exercise in 
coordinate transformation. For example, one might wish to have the incident radiation propagating along 
the (1,1,1) direction in the TF (example 14 below). Here we provide some selected examples: 

1. X = ai, y = a2, z = 33 : 8 = 0, <i> + /3 = 

2. X = ai, y = a3, z = -a2 : 8 = 0, $ + /3 = 90° 

3. X = a2, y = ai, z = -ag : 8 = 90°, f3 = 180°, $ = 

4. X = a2, y = a3, z = ai : 8 = 90°, P = 180°, $ = 90° 

5. X = a3, y = ai, z = as : 8 = 90°, 13 = -90°, $ = 
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6. X = as, y = aa, z = -ai : e = 90°, f3 = -90°, $ = -90° 

7. X = -ai, y = as, z = -as : 6 = 180°, /3 + $ = 180° 

8. i = -ai, y = aa, z = as : e = 180°, /3 + $ = 90° 

9. X = -aa, y = ai, z = as : e = 90°, /? = 0, $ = 

10. X = -as, y = ag, z = -ai : 6 = 90°, /3 = 0, $ = -90° 

11. X = -ag, y = ai, z = -aa : e = 90°, f3 = -90°, $ = 

12. X = -aa, y = aa, z = ai ; e = 90°, /? = -90°, $ = 90° 

13. i = (ai + a2)/V2, y = ag, z = (ai - a2)/V2: 9 = 45°, /3 = 180°, $ = 90° 

14. i = (ai + a2 + aa)/ V3, y = (ai - aa)/ ^2, z = (ai + aa - 2a3)/ ^6: 6 = 54.7356°, /3 = 135°, 
$ = 30°. 



17.3 Sampling in 0, $, and /3 

The present version, DDSCAT 6.1, chooses the angles /?, 6, and $ to sample the intervals (BETAMI , BETAMX), 
(THETMI , THETMX) , (PHIMIN, PHIMAX), where BETAMI, BETAMX, THETMI, THETMX, PHIMIN, 
PHIMAX are specified in ddscat . par . The prescription for choosing the angles is to: 

• uniformly sample in /3; 

• uniformly sample in $; 

• uniformly sample in cos Q. 

This prescription is appropriate for random orientation of the target, within the specified limits of /3, $, 
ande. 

Note that when DDSCAT 6.1 chooses angles it handles f3 and $ differently from 00 The range for /3 
is divided into NBETA intervals, and the midpoint of each interval is taken. Thus, if you take BETAMI=0, 
BETAMX=90, NBETA=2 you will get /3 = 22.5° and 67.5°. Similarly, if you take PHIMIN=0, PHIMAX=180, 
NPHI=2 you will get $ = 45° and 135°. 

Sampling in O is done quite differently from sampling in /3 and $. First, as akeady mentioned above, 
DDSCAT 6.1 samples uniformly in cos Q, not Q. Secondly, the sampling depends on whether NTHETA 
is even or odd. 

• If NTHETA is odd, then the values of 8 selected include the extreme values THETMI and THETMX; 
thus, THETMI=0, THETMX=90, NTHETA=3 will give you 8 = 0, 60°, 90°. 

• If NTHETA is even, then the range of cos 8 will be divided into NTHETA intervals, and the mid- 
point of each interval will be taken; thus, THETMI=0, THETMX=90, NTHETA=2 will give you 
8 = 41.41° and 75.52° [cos 8 = 0.25 and 0.75]. 

The reason for this is that if odd NTHETA is specified, than the "integration" over cos 8 is performed 
using Simpson's rule for greater accuracy. If even NTHETA is specified, then the integration over cos 8 
is performed by simply taking the average of the results for the different 8 values. 

If averaging over orientations is desired, it is recommended that the user specify an odd value of 
NTHETA so that Simpson's rule will be employed. 



18 Orientational Averaging 

DDSCAT has been constructed to facilitate the computation of orientational averages. How to go about 
this depends on the distribution of orientations which is applicable. 



This is a change from DDSCAT.4a! !. 
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18.1 Randomly-Oriented Targets 

For randomly-oriented targets, we wish to compute the orientational average of a quantity &, $): 

iQ) = -^J dp J dcosO J d$Q(/3,e,$) . (30) 

To compute such averages, all you need to do is edit the file ddscat .par so that DDSCAT knows 
what ranges of the angles (3, Q, and $ are of interest. For a randomly-oriented target with no symmetry, 
you would need to let /? run from to 360°, 6 from to 180°, and <!> from to 360°. 

For targets with symmetry, on the other hand, the ranges of /3, 8, and <!> may be reduced. First 
of all, remember that averaging over <i> is relatively "inexpensive", so when in doubt average over to 
360°; most of the computational "cost" is associated with the number of different values of which 
are used. Consider a cube, for example, with axis ai normal to one of the cube faces; for this cube [3 
need run only from to 90°, since the cube has fourfold symmetry for rotations around the axis ai. 
Furthermore, the angle 8 need run only from to 90°, since the orientation (/3,8,<E>) is indistinguishable 
from 180° - 8, 360° - 

For targets with symmetry, the user is encouraged to test the significance of /3,8,<I> on targets with 
small numbers of dipoles (say, of the order of 100 or so) but having the desired symmetry. 

It is important to remember that DDSCAT .4b handled even and odd values of NTHETA differently - 
see |8] above! For averaging over random orientations odd values of NTHETA are to be preferred, as this 
will allow use of Simpson's rule in evaluating the "integral" over cos 8. 

18.2 Nonrandomly-Oriented Targets 

Some special cases (where the target orientation distribution is uniform for rotations around the x axis = 
direction of propagation of the incident radiation), one may be able to use DDSCAT 6.1 with appropriate 
choices of input parameters. More generally, however, you will need to modify subroutine ORIENT to 
generate a hst of NBETA values of /?, NTHETA values of 8, and NPHI values of $, plus two weight- 
ing arrays WGTA ( 1-NTHETA, 1-NPHI ) and WGTB { 1-NBETA) . Here WGTA gives the weights which 
should be attached to each (8,$) orientation, and WGTB gives the weight to be attached to each /3 orien- 
tation. Thus each orientation of the target is to be weighted by the factor WGTA x WGTB. For the case of 
random orientations, DDSCAT 6.1 chooses 8 values which are uniformly spaced in cos 8, and (3 and $ 
values which are uniformly spaced, and therefore uses uniform weights 
WGTB=1./NBETA 

When NTHETA is even, DDSCAT sets 

WGTA=1./(NTHETAXNPHI) 

but when NTHETA is odd, DDSCAT uses Simpson's rule when integrating over 8 and 
WGTA= (1/3 or 4/3 or2/3)/(NTHETAxNPHI) 
Note that the program structure of DDSCAT may not be ideally suited for certain highly oriented 
cases. If, for example, the orientation is such that for a given $ value only one 8 value is possible 
(this situation might describe ice needles oriented with the long axis perpendicular to the vertical in the 
Earth's atmosphere, illuminated by the Sun at other than the zenith) then it is foolish to consider all the 
combinations of 8 and $ which the present version of DDSCAT is set up to do. We hope to improve 
this in a future version of DDSCAT. 

19 Target Generation 

DDSCAT contains routines to generate dipole arrays representing targets of various geometries, includ- 
ing spheres, ellipsoids, rectangular soUds, cylinders, hexagonal prisms, tetrahedra, two touching ellip- 
soids, and three touching ellipsoids. 
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The target type is specified by variable CSHAPE on line 9 of ddscat .par, and the target di- 
mensions (in units of interdipole spacing d) is specified by up to 6 target shape parameters (SHPARl, 
SHPAR2, SHPAR3, ...) on line 10. 

The target geometry is most conveniently described in a coordinate system attached to the target 
which we refer to as the "Target Frame" (TF), so in this section only we will let x,y,z be coordinates 
in the Target Frame. Once the target is generated, the orientation of the target in the Lab Frame is 
accomplished as described in f}T7] 

Target geometries currently supported include: 

19.1 ELLIPS = Homogeneous, isotropic ellipsoid. 

"Lengths" SHPARl, SHPAR2, SHPAR3 in the x, y, z directions in the TF. (a;/SHPARl)2+(?//SHPAR2)2+ 

(z/SHPAR3)^ = (i^/4, where d is the interdipole spacing. 

The target axes are set to ai = (1, 0, 0) and 3.2 = (0, 1, 0) in the TF. 

User must set NC0MP=1 on line 9 of ddscat .par. 

A homogeneous, isotropic sphere is obtained by setting SHPARl = SHPAR2 = SHPAR3 = diameter/d. 



19.2 ANIELL = Homogeneous, anisotropic ellipsoid. 

SHPARl, SHPAR2, SHPAR3 have same meaning as for ELLIPS. Target axes ai = (1,0, 0)anda2 = 
(0, 1,0) in the TF. It is assumed that the dielectric tensor is diagonal in the TF. User must set NC0MP=3 
and provide xx, yy, zz elements of the dielectric tensor 

19.3 CONELL = Two concentric ellipsoids 

SHPARl, SHPAR2, SHPAR3 specify the lengths along x, y, z axes (in the TF) of the outer ellipsoid; 
SHPAR4, SHPAR5, SHPAR6 are the lengths, along the x, y, z axes (in the TF) of the inner ellipsoid. 
The "core" within the inner ellipsoid is composed of isotropic material 1 ; the "mantle" between inner 
and outer ellipsoids is composed of isotropic material 2. Target axes ai = (1, 0, 0), a2 = (0, 1, 0) in TF. 
User must set NC0MP=2 and provide dielectric functions for "core" and "mantle" materials. 

19.4 CYLNDR = Homogeneous, isotropic cylinder 

Length/d=SHPARl, diameter/(i=SHPAR2, with cylinder axis = ai = (1, 0, 0) and a2 = (0, 1, 0) in the 
TF. User must set NC0MP=1. 

19.5 CYLCAP = Homogeneous, isotropic cylinder with hemispherical end-caps 

Cylinder length (not including caps!)/(i=SHPARl, diameter/rf=SHPAR2, with cylinder axis = ai = 
(1, 0, 0) and a2 = (0, 1, 0) in the TF. Note that the target length (including end-caps) will be (SHPARl 
+ SHPAR2)xd. User must setNCOMP=L 

19.6 UNICYL = Homogeneous cylinder with unixial anisotropic dielectric tensor 

SHPARl, SHPAR2 have same meaning as for CYLNDR. Cylinder axis = ai = (1,0, 0), a2 = (0, 1,0). 
It is assumed that the dielectric tensor e is diagonal in the TF, with Cyy = Czz- User must set NC0MP=2. 
Dielectric function 1 is for E jj ai (cylinder axis), dielectric function 2 is for E _L ai. 

19.7 RCTNGL = Homogeneous, isotropic, rectangular solid 

X, y, z lengths/d = SHPARl, SHPAR2, SHPAR3. Target axes ai = (1,0, 0) and a2 = (0, 1, 0) in the TF. 
User must set NC0MP=1. 
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19.8 ANIREC = Homogeneous, anisotropic, rectangular solid 

X, y, z lengths/d = SHPARl, SHPAR2, SHPAR3. Target axes ai (1,0,0) and a.2 = (0, 1,0) in the 
TF. Dielectric tensor is assumed to be diagonal in the target frame. User must set NC0MP=3 and supply 
names of three files for e as a function of wavelength or energy: first for exx, second for eyy, and third 
for e^z, as in the following sample ddscat . par file: 

' =================== Parameter file ===================' 

'**** PRELIMINARIES 

'NOTORQ'= CMT0RQ*6 (DOTORQ, NOTORQ) — either do or skip torque calculations 

'PBCGST'= CMDS0L*6 (PBCGST, PETRKP) — select solution method 

'GPFAFT'= CMETHD*6 (GPFAFT, FFTW21, CONVEX) 

' LATTDR' = CALPHA*6 (LATTDR is only supported option now) 

'NOTBIN'= CBINFLAG (ALLBIN, ORIBIN, NOTBIN) 

'NOTCDF'= CNETFLAG (ALLCDF, ORICDF, NOTCDF) 

'ANIREC' = CSHAPE*6 (FRMFIL, ELLIPS , CYLNDR, RCTNGL, HEXGON, TETRAH, UNICYL, . . . ) 

10 25 50 = shape parameters PARI, PAR2, PAR3 

3 = NCOMP = number of dielectric materials 

'TABLES'= CDIEL*6 ( TABLES , H20ICE , H20LIQ; if TABLES, then filenames follow...) 

' /u/draine/DDA/diel/m2 . 00_0 . 10' = name of file containing dielectric function 

' /u/draine/DDA/diel/ml . 50_0 . 00' 

' /u/draine/DDA/diel/ml . 50_0 . 00' 

'**** CONJUGATE GRADIENT DEFINITIONS ****' 

= INIT (TO BEGIN WITH | X0> = 0) 

l.OOe-5 = ERR = MAX ALLOWED (NORM OF | G>=AC I E>-ACA | X> ) / (NORM OF AC|E>) 
' -k-k-k-k Angular resolution for calculation of <cos>, etc. + + + + ' 
2.0 = ETASCA (number of angles is proportional to [ (2+x) /ETASCA] "2 ) 

'-k-k-k-k Wavelengths (micron) -k-k-k-k' 

500. 500. 1 ' INV = wavelengths ( first , last , how many , how=LIN, INV, LOG) 
' -k-k-k-k Effective Radii (micron) ' 

30.60 30.60 1 'LIN' = eff. radii (first, last, how many, how=LIN, INV, LOG) 
' kkkk Define Incident Polarizations ****' 

(0,0) (l.,0.) (0.,0.) = Polarization state eOl (k along x axis) 

2 = lORTH (=1 to do only pol . state eOl; =2 to also do orth. pol . state) 

1 = IWRKSC (=0 to suppress, =1 to write ".sea" file for each target orient. 
'-k-k-k-k Prescribe Target Rotations ****' 

0. 0. 1 = BETAMI, BETAMX, NBETA (beta=rotation around al) 
0. 90. 3 = THETMI, THETMX, NTHETA (theta=angle between al and k) 
0. 0. 1 = PHIMIN, PHIMAX, NPHI (phi=rotation angle of al around k) 
'**** Specify first IWAV, TRAD, lORI (normally 0) ****' 
= first IWAV, first IRAD, first TORI (0 to begin fresh) 
'-k-k-k-k Select Elements of S_ij Matrix to Print ****' 

6 = NSMELTS = number of elements of S_ij to print (not more than 9) 

11 12 21 22 31 41 = indices ij of elements to print 
'-k-k-k-k Specify Scattered Directions ****' 

0. 0. 180. 30 = phi, thetan_min, thetan_max, dtheta (in degrees) for plane A 
90. 0. 180. 30 = phi, ... for plane B 

19.9 HEXGON = Homogeneous, isotropic hexagonal prism 

SHPARl = (Length of prism)/c? = (distance between hexagonal faces)/(i, SHPAR2 = (distance between 
opposite vertices of one hexagonal face)/o? = 2 x hexagon side/d. Note that Target axis ai = (1, 0, 0) in 
the TF is along the prism axis, and target axis 3.2 = (0, 1 , 0) in the TF is normal to one of the rectangular 
faces of the hexagonal prism. User must set NC0MP=1. 

19.10 TETRAH = Homogeneous, isotropic tetrahedron 

SHPARl=length/(i of one edge. Orientation: one face parallel to y,z plane (in the TF), opposite "vertex" 
is in +x direction, and one edge is parallel to z axis (in the TF). Target axes ai = (1, 0, 0) [emerging 
from one vertex] and a2 = (0, 1, 0) [emerging from an edge] in the TF. User must set NC0MP=1. 
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19.11 TWOSPH = Two touching homogeneous, isotropic spheroids, with distinct 
compositions 

First spheroid has length SHPARl along symmetry axis, diameter SHPAR2 perpendicular to symmetry 
axis. Second spheroid has length SHPAR3 along symmetry axis, diameter SHPAR4 perpendicular to 
symmetry axis. Contact point is on line connecting centroids. Line connecting centroids is in x direction. 
Symmetry axis of first spheroid is in y direction. Symmetry axis of second spheroid is in direction 
y cos(SHPAR5) + z sin(SHPAR5), where y and z are basis vectors in TF, and SHPAR5 is in degrees. If 
SHPAR6=0., then target axes ai = (1, 0, 0), a2 = (0, 1, 0). If SHPAR6=1., then axes ai and a2 are set 
to principal axes with largest and 2nd largest moments of inertia assuming spheroids to be of uniform 
density. User must set NC0MP=2 and provide dielectric function files for each spheroid. 

19.12 TWOELL = Two touching, homogeneous, isotropic ellipsoids, with distinct 
compositions 

SHPARl, SHPAR2, SHPAR3=x-length/d, y-length/d, z-length/d of one ellipsoid. The two ellipsoids 
have identical shape, size, and orientation, but distinct dielectric functions. The line connecting ellipsoid 
centers is along the a--axis in the TF. Target axes ai = (1, 0, 0) [along line connecting ellipsoids] and 
3.2 = (0, 1, 0). User must set NC0MP=2 and provide dielectric function file names for both ellipsoids. 
Ellipsoids are in order of increasing x: first dielectric function is for ellipsoid with center at negative x, 
second dielectric function for ellipsoid with center at positive x. 

19.13 TWOAEL = Two touching, homogeneous, anisotropic ellipsoids, with dis- 
tinct compositions 

Geometry as for TWOELL; SHPARl, SHPAR2, SHPAR3 have same meanings as for TWOELL. Target 
axes ai ~ (1, 0, 0) and 3.2 = (0, 1, 0) in the TF. It is assumed that (for both ellipsoids) the dielectric 
tensor is diagonal in the TF. User must set NC0MP=6 and provide xx, yy, zz components of dielectric 
tensor for first ellipsoid, and xx, yy, zz components of dielectric tensor for second ellipsoid (ellipsoids 
are in order of increasing x). 

19.14 THRELL = Three touching homogeneous, isotropic ellipsoids of equal size 
and orientation, but distinct compositions 

SHPARl, SHPAR2, SHPAR3 have same meaning as for TWOELL. Line connecting ellipsoid centers is 
parallel to x axis. Target axis ai — (1, 0, 0) along line of ellipsoid centers, 3.2 = (0, 1, 0). User must set 
NC0MP=3 and provide (isotropic) dielectric functions for first, second, and third ellipsoid. 

19.15 THRAEL = Three touching homogeneous, anisotropic ellipsoids with same 
size and orientation but distinct dielectric tensors 

SHPARl, SHPAR2, SHPAR3 have same meanings as for THRELL. Target axis ai = (1,0,0) along 
line of ellipsoid centers, 3.2 = (0, 1, 0). It is assumed that dielectric tensors are all diagonal in the TF. 
User must set NC0MP=9 and provide xx, yy, zz elements of dielectric tensor for first ellipsoid, xx, yy, 
zz elements for second ellipsoid, and xx, yy, zz elements for third elhpsoid (elhpsoids are in order of 
increasing x). 

19.16 BLOCKS = Homogeneous target constructed from cubic "blocks" 

Number and location of blocks are specified in separate file blocks . par with following structure: 
one line of comments (may be blank) 
PRIN (= or 1 - see below) 
N (= number of blocks) 
B (= width/d of one block) 



TARGET GENERATION 



X y z (= position of 1st block in units of Bd) 
X y z {= position of 2nd block) 

X y z (= position of Nth block) 
If PRIN=0, then ai = (1, 0, 0), a2 = (0, 1, 0). If PRIN=1, then ai and a2 are set to principal axes with 
largest and second largest moments of inertia, assuming target to be of uniform density. User must set 

NC0MP = 1. 

19.17 DW1996 = 13 block target used by Draine & Weingartner (1996). 

Single, isotropic material. Target geometry was used in study by Draine & Weingartner (1996) of radia- 
tive torques on irregular grains, ai and a2 are principal axes with largest and second-largest moments of 
inertia. User must set NC0MP = 1. Target size is controlled by shape parameter SHPAR { 1 ) = width of 
one block in lattice units. 

19.18 NSPHER = Multisphere target consisting of the union of spheres of 
single isotropic material 

Spheres may overlap if desired. The relative locations and sizes of these spheres are defined in an 
external file, whose name (enclosed in quotes) is passed through ddscat . par. The length of the file 
name should not exceed 13 characters. The pertinent line in ddscat. par should read 
SHPAR (1) SHPAR (2) '_^/e«ame' (quotes must be used) 

where SHPAR ( 1 ) = target diameter in x direction (in Target Frame) in units of d 
SHPAR (2) = to have ai = (1,0, 0), 02 = (0, 1,0) in Target Frame. 

SHPAR ( 2 ) = 1 to use principal axes of moment of inertia tensor for oi (largest /) and 02 (intermediate 
/). 

filename is the name of the file specifying the locations and relative sizes of the spheres. 
The overall size of the multisphere target (in terms of numbers of dipoles) is determined by parameter 
SHPAR ( 1 ) , which is the extent of the multisphere target in the x-direction, in units of the lattice spacing 
d. The file filename ' should have the following structure: 

N (= number of spheres) 

one line of comments (may be blank) 

xi yi z\ a\ (arb. units) 

X2 1/2 Z2 02 (arb. units) 

'^N yN zn un (arb. units) 
where the coordinates of the center, and Uj is the radius of sphere j. 

Note that xj, yj, Zj, aj (j = 1, N) establish only the shape of the A^— sphere target. For instance, 
a target consisting of two touching spheres with the line between centers parallel to the x axis could 
equally well be described by lines 3 and 4 being 

0000.5 

1 000.5 

or 

000 1 
200 1 

The actual size (in physical units) is set by the value of Ocff specified in ddscat . par, where, as always, 
Coff = (3y/47r)^/'^, where V is the volume of material in the target. 
User must set NC0MP=1. 

19.19 PRISM3 = Triangular prism of homogeneous, isotropic material 

SHPARl, SHPAR2, SHPAR3, SHPAR4 = a/rf, &/a, c/a, i/a 

The triangular cross section has sides of width a, b, c. L is the length of the prism, d is the lattice 
spacing. The triangular cross-section has interior angles a, (3, 7 (opposite sides a, 6, c) given by cos a — 
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(6^ + - a'^)/2bc, cos/3 = (a^ + - b'^)/2ac, cos 7 = (a^ + 6^ - c^)/2ab. In the Target Frame, 
the prism axis is in the x direction, the normal to the rectangular face of width a is (0,1,0), the normal to 
the rectangular face of width b is (0, — cos 7, sin 7), and the normal to the rectangular face of width c is 

(0, — cos 7, — sin 7). 



19.20 LYRSLB = Multilayer rectangular slab 

Multilayer rectangular slab with overall x, y, z lengths = SHPARlxd, ay = SHPAR2x(i, = 
SHPARSxd, with the boundary between layers 1 and 2 occuring at xi = SHPAR4 x a.j., boundary 
between layers 2 and 3 occuring at X2 = SHPAR5 x a.^. [if SHPAR5> 0], and boundary between lay- 
ers 3 and 4 occuring at .T3 = SHPAR6 x [if SHPAR6> 0], To create a simple bilayer slab, just set 
SHPAR5=SHPAR6=0. User must set NCOMP=2,3, or 4 and provide dielectric function files for each of 
the two layers. 



19.21 FRMFIL = Target defined by list of dipole locations and "compositions" 
obtained from a file 

This option causes DDSCAT to read the target geometry information from a file shape . dat instead 
of automatically generating one of the geometries listed above. The shape . dat file is read by routine 
REASHP (file reashp . f ). The user can customize REASHP as needed to conform to the manner in 
which the target geometry is stored in file shape . dat. However, as supphed, REASHP expects the file 
shape . dat to have the following structure: 

• one line containing a description; the first 67 characters will be read and printed in various output 
statements. 

• N = number of dipoles in target 

• o^ix aiy aiz = x,y,z components of ai 

• 0i2x (i2y o,2z = x,y,z Components of a2 

• (line containing comments) 

• dummy IXYZ {1, 1) IXYZ(1,2) IXYZ(1,3) IC0MP(1,1) IC0MP(1,2) IC0MP(1,3) 

• dummy IXYZ (2, 1) IXYZ (2, 2) IXYZ (2, 3) IC0MP(2,1) ICOMP(2,2) ICOMP(2,3) 

• dummy IXYZ (3, 1) IXYZ (3, 2) IXYZ (3, 3) ICOMP (3, 1) ICOMP(3,2) ICOMP(3,3) 

• ... 

• dummy IXYZ {J, 1) IXYZ (J, 2) IXYZ (J, 3) ICOMP (J, 1) ICOMP (J, 2) ICOMP (J, 3) 

• ... 

• dtimmy IXYZ (N, 1 ) IXYZ (N, 2) IXYZ (N, 3) ICOMP (N,l) ICOMP (N, 2) ICOMP (N, 3) 



19.22 Modifying Existing Routines or Writing New Ones 

The user should be able to easily modify these routines, or write new routines, to generate targets with 
other geometries. The user should first examine the routine target . f and modify it to call any new 
target generation routines desired. Alternatively, targets may be generated separately, and the target 
description (locations of dipoles and "composition" corresponding to x,y,z dielectric properties at each 
dipole site) read in from a file by invoking the option FRMFIL in ddscat . f . 

It is often desirable to be able to run the target generation routines without running the entire 
DDSCAT code. We have therefore provided a program CALLTARGET which allows the user to generate 
targets interactively; to create this executable just tvpJ*! 

" Non-Unix sites: The source code for CALLTARGET is in the file MISC . FOR. You must compile and link this to ERRMSG, 
GETSET, REASHP, TAR2EL, TAR2SP, TAR3EL, TARBLOCKS, TARCEL, TARCYL, TARELL, TARGET, TARHEX, TARREC, 
TARTET, and WRIMSG. The source code for ERRMSG is in SRC2 . FOR, GETSET is in SRC5 . FOR, and the rest are in SRC8 . FOR 
and SRC 9. FOR. 
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make calltarget . 

The program calltarget is to be run interactively; the prompts are self-explanatory. You may need 
to edit the code to change the device number IDVOUT as for DDSCAT (see ^6.2.1l above). 

After running, calltarget will leave behind an ASCII file target . out which is a list of the 
occupied lattice sites in the last target generated. The format of target . out is the same as the format 
of the shape . dat files read if option FRMFIL is used (see above). Therefore you can simply 

mv target. out shape.dat 
and then use DDSCAT with the option FRMFIL in order to input a target shape generated by CALLTARGET. 

20 Scattering Directions 

DDSCAT calculates scattering in selected directions, and elements of the scattering matrix are reported 
in the output files vxxryykzzz . sea . The scattering direction is specified through angles Og and (pg (not 
to be confused with the angles Q and $ which specify the orientation of the target relative to the incident 
radiation!). 

The scattering angle 6s is simply the angle between the incident beam (along direction x) and the 
scattered beam {6s — for forward scattering, 9s = 180° forbackscattering). 

The scattering angle (ps specifies the orientation of the "scattering plane" relative to the ic, y plane 
in the Lab Frame. When <j>s = the scattering plane is assumed to coincide with the ic, y plane. When 
(ps = 90° the scattering plane is assumed to coincide with the x, z plane. Within the scattering plane the 
scattering directions are specified by < 6*^ < 180°. 

Scattering directions for which the scattering properties are to be calculated are set in the parameter 
file ddscat . par by specifying one or more scattering planes (determined by the value of </>«) and for 
each scattering plane, the number and range of 6s values. The only limitation is that the number of 
scattering directions not exceed the parameter MXSCA in DDSCAT . f (in the code as distributed it is set 
to MXSCA=1000). 

21 Incident Polarization State 

Recall that the "Lab Frame" is defined such that the incident radiation is propagating along the x axis. 
DDSCAT 6.1 allows the user to specify a general elliptical polarization for the incident radiation, by 
specifying the (complex) polarization vector Gqi. The orthonormal polarization state eo2 = x x e^-^ is 
generated automatically if ddscat .par specifies I0RTH=2. 

For incident linear polarization, one can simply set Gqi = y by specifying (0,0) (1,0) (0,0) 
in ddscat . par; then eo2 = z For polarization mode Gqi to correspond to right-handed circular po- 
larization, set Gqi = (y + iz)/-y/2 by specifying (0,0) (1,0) (0,1) in ddscat . par (DDSCAT 
6.1 automatically takes care of the normalization of Gqi); then eo2 = (— ty + z)/-y/2, corresponding to 
left-handed circular polarization. 

22 Averaging over Scattering Directions: g{l) = (cos 9s), etc. 
22.1 Angular Averaging 

An example of scattering by a nonspherical target is shown in Fig. [8] showing the scattering for a tilted 
cube. Results are shown for four different scattering planes. 

DDSCAT automatically carries out numerical integration of various scattering properties, including 

• (cos 6*5 ); 

• (cos^ 6's); 

• g = (cos 6s)x+ (sin 9s cos <j>s)y + (sin 9s sin (f>s)z (see ^ 031 ); 

• Qr, provided option DOTORQ is specified (see ^T5[ . 
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Figure 8: Scattered intensity for incident unpolarized light for a cube with m = 1.5 + O.OOli and x — 27raoft/A = 6.5 (i.e., 
d/X — 1.6676, where d is the length of a side). The cube is tilted with respect to the incident radiation, with O = 30°, and rotated 
by /3 = 15° around its axis to break reflection symmetry. The Mueller matrix element Sn is shown for 4 different scattering planes. 
The strong forward scattering lobe is evident. It is also seen that the scattered intensity is a strong function of scattering angle (ps as 
well as 9s - at a given value of 63 (e.g., 9s = 120°), the scattered intensity can vary by orders of magnitude as (j>s changes by 90°. 



The angular averages are accomplished by evaluating the scattered intensity for selected scattering di- 
rections {6s, <f>s), and taking the appropriately weighted sum. Suppose that we have Ng different values 
of 9 s, 

d = dj, j^l,...,Ne , (31) 
and for each value of 9j, N^{j) different values of 



t>j,k, k = l,...,N^{j) 



(32) 



For a given j, the values of (pj^^ are assumed to be uniformly spaced: 4>j_k+i ~ 4'j,k = 27r/iV0(j). The 
angular average of a quantity f{Os, 4>s) is approximated by 



1 

~ in 



sin OsdOs 



d(i)sf{9s,4's 



, No N4,{3) 

[zE E /(^.'^..^)". 



j=i k=i 



Ng,k 



[cos(6ij„i) -cos(6ij+i)] , j=2,...,Ne-l 



27r 



27r 



cos(0i) + cos(6'2) 



For a sufficiently large number of scattering directions 



1 



COs(6'Ar^_l) + COs(6ljvJ ' 

2 



Ng 

sea -^^N^ij) 
J = l 



(33) 

(34) 
(35) 
(36) 

(37) 



the sum ( |33] l approaches the desired integral, but the calculations can be a significant cpu-time burden, 
so efficiency is an important consideration. 
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22.2 Selection of Scattering Angles Og , (/^s 

Beginning with DDSCAT 6.1, an improved approach is taken to evaluation of the angular average. Since 
targets with large values of x = 27racff / A in general have strong forward scattering, it is important to 
obtain good sampling of the forward scattering direction. To implement a preferential sampling of the 
forward scattering directions, the scattering directions {6, (f>) are chosen so that 9 values correspond 
to equal intervals in a monotonically increasing function s{9). The function s{d) is chosen to have a 
negative second derivative cPs/dO'^ < so that the density of scattering directions will be higher for 
small values of 0. DDSCAT 6.1 takes 

This provides increased resolution (i.e., increased ds/dO) for 6 < 6q. We want this for the forward 
scattering lobe, so we take 

Oo = . (39) 

1 + a; 

The 9 values run from 9 = to 9 = tt, corresponding to uniform increments 



If we now require that 



A. = ^U + — ^) . (40) 



r . .1 As 7r/2 
max[A0]«— — =Vir^ , (41) 



we determine the number of values of 9: 

2(3 + x) [l + l/in + 9o)] 

V [i + V(- + ^o)^] ■ ^ ^ 

Thus for small values of x, max[A6'] 30°7^ and for x ^ 1, max[A0] 90°ri/x. For a sphere, 
minima in the scattering pattern are separated by ^ 180°/a:;, so 77 = 1 would be expected to marginally 
resolve structure in the scattering function. Smaller values of 77 will obviously lead to improved sampUng 
of the scattering function, and more accurate angular averages. 

The scattering angles 9j used for the angular averaging are then given by 



{s, - 1 - 9o) + [{1 + 9o - s,)^ + i9o 



1/2 



2 

where 



(43) 



Ne-1 



1 + ^ 



j = l,...,Ne . (44) 



For each 9j, we must choose values of </>. For 9i = Q and 9^^ = tt only a single value of (p is needed 
(the scattering is independent of (p in these two directions). For < 9j < tt we use 

= max{3,nint [47rsin(6lj)/(e'j+i - 9j-i)]} , (45) 

where nint = nearest integer. This provides sampling in (f> consistent with the sampling in 9. 



22.3 Accuracy of Angular Averaging as a Function of r] 

Figure |9] shows the absolute errors in {cos 9) and (cos^ 9) calculated for a sphere with refractive index 
m = 1.33 + O.Olz using the above prescription for choosing scattering angles. The error is shown as a 
function of scattering parameter x. We see that accuracies of order 0.01 are attained with rj = 1, and that 
the above prescription provides an accuracy which is approximately independent of x. We recommend 
using values of 77 < 1 unless accuracy in the angular averages is not important. 
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Figure 9: Errors in (cos 6) and (cos'^ 9) calculated for z N — 17904 dipole pseudosphere with m — 1.33 + O.Oli, as functions 
of a:: = 2na/X. Results are shown for different values of the parameter 77. 
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Figure 10: Same as Fig.[9l but f or a Af = 32768 dipole cube with m = 1.5 + O.OOli. 
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23 Mueller Matrix for Scattering in Selected Directions 
23.1 Two Orthogonal Incident Polarizations (I0RTH=2) 

DDSCAT 6.1 internally computes the scattering properties of the dipole array in terms of a complex 
scattering matrix fmii^s, <l>s) (Draine 1988), where index 1 — 1,2 denotes the incident polarization state, 
m = 1,2 denotes the scattered polarization state, and 6s,<f>s specify the scattering direction. Normally 
DDSCAT is used with I0RTH=2 in ddscat .par, so that the scattering problem will be solved for 
both incident polarization states (1 = 1 and 2); in this subsection it will be assumed that this is the case. 

Incident polarization states / = 1, 2 correspond to polarization states eoi, eo2; recall that polarization 
state eoi is user-specified, and eo2 = x x e^^. Scattered polarization state m = 1 corresponds to 
linear polarization of the scattered wave parallel to the scattering plane (ei = ep^ = 9s) and m = 2 
corresponds to linear polarization perpendicular to the scattering plane (in the +(j)s direction). The 
scattering matrix /,„/ was defined (Draine 1988) so that the scattered electric field Es is related to the 
incident electric field (0) at the origin (where the target is assumed to be located) by 



E, • , 



exp(ik • r) 



hi 
hi 



/l2 
/22 



E,(0) - eoi 

E,(0) -602 



(46) 



The 2x2 complex scattering amplitude matrix (with elements 5*1, 5*2, 5*3, and S4) is defined so that (see 
Bohren & Huffman 1983) 



E, • 

-E, • i 



cxp(ik • r) 
—ikr 



S2 S3 
S4. Si 



E.(0) 
E,(0) ■ 



e.-i|| 



(47) 



where eq , e^j^ are (real) unit vectors for incident polarization parallel and perpendicular to the scattering 



plane (with the customary definition of e, j 
From ( |46|47| i we may write 



= e,|| X x). 



^2 
^4 



^3 
^1 



E,(0) 
E.(0) . 



fll fl2 
^hl ^./22 



E.(0) 
E.(0) 



-02 



(48) 



Let 



a = 


^01 


•y , 


(49) 


b = 


^01 


•z , 


(50) 


c = 


^02 


•y , 


(51) 


d = 


^02 


• z . 


(52) 



Note that since eoi, eo2 could be complex (i.e., elliptical polarization), the quantities a, 6, c, d are com- 
plex. Then 

eSn / \ c d j \ i 



(53) 



and eq. ( |48T l can be written 



'S'2 <S'3 

tS'4 5*1 



E.(0)-e,|| 
E,(0)-e,i 



-hi 

hi 



-/12 

/22 



E.(0)-y 
E,(0) -z 



(54) 



The incident polarization states 6^11 and e^j^ are related to y, z by 



substituting dSSt into ( |54] | we obtain 



'S'2 <5'3 

Si Si 



E,(0) -e,,! 
E.(0) -e^i 



cos ( 
sin ( 



-hi 
hi 



-/12 

/22 



sm q>. 

~ COS( 



a h 
c d 



cos ( 
sinf 



sm (p; 

- COS( 



(55) 



E,(0) • e,|| \ 
E,(0) •e.i ) 
(56) 
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Eq. ( |56] | must be true for all Ei(0), so we obtain an expression for the complex scattering amplitude 
matrix in terms of the f„ii : 

f ? ? ) - ^ f V" )(" """t "'"t ) ■ (57) 

\^ ^4 5^1 / V /21 h2 ) \ c d J \ s\ncj)s -coscj)s J 

This provides the 4 equations used in subroutine GETMUELLER to compute the scattering amplitude 
matrix elements: 

51 = -i [f 2i{b cos 4>s - a sin 4>s) + f 22{d cos 4>s - c sin 4>s)] , (58) 

52 = -i [/ii(acos0s + 6sin0s) + /i2(ccos0s + dsin^s)] , (59) 

53 = i [/ii(6cos0s - asin^s) + /i2(dcos0s - csin^s)] , (60) 
<S'4 = i [f 2iia cos (t)s + b sin (l)s) + f 22{c cos (j)s + d sink's)] . (61) 



23.2 Stokes Parameters 

It is both convenient and customary to characterize both incident and scattered radiation by 4 "Stokes 
parameters" - the elements of the "Stokes vector". There are different conventions in the literature; we 
adhere to the definitions of the Stokes vector {I,Q,U,V) adopted in the excellent treatise by Bohren & 
Huffman (1983), to which the reader is referred for further detail. Here are some examples of Stokes 
vectors (/, Q, C/, V) = (1, Q//, U/I, V/I)I: 

• (1,0, 0,0)/: unpolarized light (with intensity /); 

• (1,1, 0, 0)/ : 100% linearly polarized with E parallel to the scattering plane; 

• (1,-1,0, 0)/ : 100% linearly polarized with E perpendicular to the scattering plane; 

• (1, 0, 1, 0)/ : 100% linearly polarized with E at +45° relative to the scattering plane; 

• (1,0,-1, 0)/ : 100% linearly polarized with E at -45° relative to the scattering plane; 

• (1,0, 0, 1)/ : 100% right circular polarization {i.e., negative helicity); 

• (1, 0, 0, —1)/ : 100% left circular polarization (i.e., positive helicity). 



23.3 Relation Between Stokes Parameters of Incident and Scattered Radiation: 
The Mueller Matrix 

It is convenient to describe the scattering properties in terms of the 4x4 Mueller matrix relating the 
Stokes parameters {li, Qi,Ui,Vi) and {Is, Qs, Us, Vs) of the incident and scattered radiation: 







/ ^11 


Sl2 


Sl3 


Si4 
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I \ 


Qs 


1 


'S'21 


S22 
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S2i 




Q^ 
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Ur 


V Vs J 
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S'42 


'S'43 


S'44 







(62) 



Once the amplitude scattering matrix elements are obtained, the Mueller matrix elements can be com- 
puted (Bohren & Huffman 1983): 



Sii 


- (|5lP + |52p + |53p + 


\S4\') /2 


S12 


- {\S2?~\Sl? + \S,\'^ 


\S3\') /2 


Sl3 


= Re{S2S; + SiSl) , 




SlA 


= lmiS2S*s~SiSl) , 




S21 


= (1^21' -|5i|' + 1^31'- 


\S,f) /2 


S22 


= (|5lp + |52p-|53p- 


\S,\') /2 


S23 


= Re (S'25'3 — S'lS'l) , 
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S24 


= lm{S2S*3 


+ ^1^4*) 


S31 


= ReiS2Sl 


+ Sis;) 


S32 


= RciS2Sl 


~ S1S3) 




= Re{SiS; 


+ S3St) 


S34 


= Im(S'2S'i 


+ 848^) 


S41 


= Im {S4S2 


+ Si 3*3) 


S42 


= Im {S4S2 


— SiS^) 


843 


= Im(5i52* 


— 5*351) 


S44 




— S'354) 



These matrix elements are computed in DDSCAT and passed to subroutine WRITESCA which handles 
output of scattering properties. Although the Muller matrix has 16 elements, only 9 are independent. 

The user can select up to 9 distinct Muller matrix elements to be printed out in the output files 
wxxryykzzz . sea and vixxryyori . avg; this choice is made by providing a list of indices in ddscat . pa 
(see AppendixlAb. 

If the user does not provide a list of elements, WRITESCA will provide a "default" set of 6 selected 
elements: Su, S21, S31, S41 (these 4 elements describe the intensity and polarization state for scattering 
of unpolarized incident radiation), S12, and 5i3. 

In addition, WRITESCA writes out the linear polarization P of the scattered light for incident unpo- 
larized light (Bohren & Huffman 1983, p. 382) 

P = -^ . (64) 



23.4 Polarization Properties of the Scattered Radiation 



The scattered radiation is fully characterized by its Stokes vector {Is,Qs,Us,Vs). As discussed in 
Bohren & Huffman (eq. 2.87), one can determine the linear polarization of the Stokes vector by op- 
erating on it by the Mueller matrix of an ideal linear polarizer: 



/ 1 



Spol 



cos 2^ sin 2^ \ 

cos 2^ 008^2^ cos 2^ sin 2^ 

sin 2^ sin 2^ cos 2^ sin^ 2^ 

/ 



(65) 



where ^ is the angle between the unit vector 6s parallel to the scattering plane ("SP") and the "transmis- 
sion" axis of the linear polarizer Therefore the intensity of light polarized parallel to the scattering plane 
is obtained by taking ^ = and operating on {Is, Qs, Us, Vs) to obtain 



/(E, II SP) 



;{Is + Qs) 



(66) 



2fc2r2 



[(5ll + 52l)/. + (^12 + 522)Q, + (5i3 + S23)U., + {Si4 + 824^] . (67) 



Similarly, the intensity of light polarized perpendicular to the scattering plane is obtained by taking 



/(E, ±SP) = l{Is-Qs) 



(68) 



2fc2r2 



[(5ll - 52l)/. + (5i2 - 522)g. + (5i3 - S23)U, + (5l4 - S24)V^ . (69) 
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23.5 Relation Between Mueller Matrix and Scattering Cross Sections 

Differential scattering cross sections can be obtained directly from the Mueller matrix elements by noting 
that 

Let SP be the "scattering plane": the plane containing the incident and scattered directions of propaga- 
tion. Here we consider some special cases: 

• Incident light unpolarized: Stokes vector = /(I, 0, 0, 0): 

- cross section for scattering with polarization Es || SP: 

- cross section for scattering with polarization -L SP: 

+ (72, 

- total intensity of scattered light: 



(l^lP + \S2\' + |^3P + 1^41') - T^^ll (73) 



• Incident light polarized with E^ || SP: Stokes vector Si = /(1, 1, 0, 0): 
- cross section for scattering with polarization Eg || to SP: 



-rfg^= p 1^21' = ^ (511+^12+521 +^22) (74) 



- cross section for scattering with polarization ± to SP: 

dCana 1 , ^ , o 1 



= pl*^*!^ = ^ (5ii +5i2 - 521 - 522) (75) 



- total scattering cross section: 



"^^"^"^ ^ {\S2? + 154^ = ^ (5ii + 5i2) (76) 



F VI -I ' I ^1 ; p 

Incident light polarized with Ei + SP: Stokes vector Si — /(I, —1, 0, 0): 
- cross section for scattering with polarization Eg || to SP: 



- cross section for scattering with polarization Eg + SP: 



2 |53p — ■;rr2 (5ii ^ 5i2 + 521 — S22) {11) 



dVL - pi-^il^ - ^(-^ii ^'^12 -521 +522) (78) 

- total scattering cross section: 

^ = ^(|5ip + |53n^^(5ii-5i2) (79) 

Incident light linearly polarized at angle 7 to SP: Stokes vector Si = /(I, cos 27, sin 27, 0): 

- cross section for scattering with polarization Eg || to SP: 



dn 2k^ 



[(5ii + 52i) + (5i2 + 522) cos 27 + (5i3 + 523) sin 27] (80) 
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- cross section for scattering with polarization _L SP: 

dCsca 1 



dn 2fc2 

- total scattering cross section: 



[(^11 - ^2i) + (5i2 - ^22)cos27+ (5i3 - ^23) Sin 27] (81) 



dn fc2 



[Sii + S12 cos 27 + 5i3 sin 27] (82) 



23.6 One Incident Polarization State Only (I ORTH= 1) 

In some cases it may be desirable to limit the calculations to a single incident polarization state - for 
example, when each solution is very time-consuming, and the target is known to have some symmetry 
so that solving for a single incident polarization state may be sufficient for the required purpose. In this 
case, set I0RTH=1 in ddscat .par. 

When I0RTH=1, only /n and /21 are available; hence, DDSCAT cannot automatically generate 
the Mueller matrix elements. In this case, the output routine WRITE SCA writes out the quantities |/ii p, 
1/21 p, Re(/ii/2i), and Im(/ii/|]^) for each of the scattering directions. 

The differential scattering cross section for scattering with polarization Eg \\ and _L to the scattering 
plane are 

'^^"^''^ Al/nl' (83) 



dn y^ ii fc2' 

Note, however, that if IPHI is greater than 1, DDSCAT will automatically set I0RTH=2 even if 
ddscat .par specified I0RTH=1: this is because when more than one value of the target orientation 
angle $ is required, there is no additional "cost" to solve the scattering problem for the second incident 
polarization state, since when solutions are available for two orthogonal states for some particular target 
orientation, the solution may be obtained for another target orientation differing only in the value of 
$ by appropriate linear combinations of these solutions. Hence we may as well solve the "complete" 
scattering problem so that we can compute the complete Mueller matrix. 



24 Using MPI (Message Passing Interface) 

Many (probably most) users of DDSCAT will wish to calculate scattering and absorption by a given 
target for more than one orientation relative to the incident radiation (e.g., when calculating orientational 
averages). DDSCAT can automatically carry out the calculations over different target orientations. Nor- 
mally, these calculations will be carried out sequentially. However, on a multiprocessor system [either 
a symmetric multiprocessor (SMP) system, or a cluster of networked computers] with MP I (Message 
Passing Interface) software installed, DDSCAT 6.1 allows the user to carry out parallel calculations of 
scattering by different target orientations, with the results gathered together and with orientational aver- 
ages calculated just as they are when sequential calculations are carried out. Note that the MPI-capable 
executable can also be used for ordinary serial calculations using a single cpu. 

More than one implementation of MPI exists. MPI support within DDSCAT 6.1 is compliant with 
MPI- 1.2 and MPI-2 standard^ and should be usable under any implementation of MPI that is compat- 
ible with those standards. 

At Princeton University Department of Astrophysical Sciences we are using MPIClQ a publicly- 
available implementation of MP I . 

http : / /www . mpi-f orum. org/| 
'' http : / /www . mcs . ani . gov/mpi/mpich| 
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24.1 Source Code for the MPI-compatible version of DDSCAT 

In order to use DDSCAT with MPI, the user must first make sure that the correct version of the code is 
being used. To prepare the MPI-capable version: 

1. In DDSCAT . f : 

• Make sure the Hne 

INCLUDE 'mpif.h' 
is not commented out. 

• Make sure the Hne 

C INTEGER MPI_COMM_WORLD 

is commented out. 

2. Inmpisubs.f, SUBROUTINE COLSUM : 

• Make sure the Hne 

INCLUDE 'mpif.h' 
is not commented out. 

• Make sure the Hne 

C INTEGER MPI_COMM_WORLD, MPI_SUM, MPI_REAL, MPI.COMPLEX 

is commented out. 

3. Inmpisubs.f, SUBROUTINE SHAREl : 

• Make sure the Hne 

INCLUDE 'mpif.h' 
is not commented out. 

• Make sure the lines 

C INTEGER MPI_CHARACTER,MPI_COMM_WORLD,MPI_INTEGER, 

C & MPI_INTEGER2,MPI_REAL,MPI_C0MPLEX 

are commented out. 

4. Inmpisubs.f, SUBROUTINE SHARE2 : 

• Make sure the Hne 

INCLUDE 'mpif.h' 
is not commented out. 

• Make sure the Hne 

C INTEGER MPI_COMM_WORLD, MPI_REAL, MPI_COMPLEX 

is commented out. 

24.2 Compiling and Linking the MPI-compatible version of DDSCAT 

In the Makefile, enable the following definitions: 

MPIsrc = mpi_subs . f 

MPIobj = mpi_subs.o 

At Princeton University Dept. of Astrophysical Sciences we are currently using mpich-1.2. 
mpich provides a "compiler wrapper" mpif77 that is supposed to automatically set the variables 
MP I CH_F 7 7 and MP I CH_F 7 7 L INKER to values appropriate for use of the pgf 7 7 compiler from PG^^ 
which is installed on our beowulf cluster. Other compilers, e.g., the Intel f77 compiler, can also be 
employed. 

24.3 Code Execution Under MPI 

Local installations of MP I will vary - you should consult with someone familiar with way MP I is in- 
stalled and used on your system. 

^ Ihttp : 77www-unix .mcs . anl . gov/mpi/mpich/| 
^' [http : / /www . pgroup . com| 
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At Princeton University Dept. of Astrophysical Sciences we use PBS (Portable Batch Systemic to 
schedule jobs. MP I jobs are submitted using PBS by first creating a shell script such as the following 
example file pbs . submit: 

# ! /bin/bash 

#PBS -1 nodes=2 :ppn=l 

#PBS -1 mem=12 00MB, pmem=300MB 

#PBS -m bea 

#PBS -j oe 

Cd $PBS_0_WORKDIR 

/usr/local/bin/mpiexec ddscat 

The lines beginning with #PBS -1 specify the required resources: 
#PBS -1 nodes=2 : ppn=l specifies that 2 nodes are to be used, with 1 processor per node. 
#PBS -1 mem=l 2 0MB, pmem=3 0MB specifies that the total memory required (mem) is 1200MB, 
and the maximum physical memory used by any single process (pmem) is 300MB. The actual definition 
of mem is not clear, but in practice it seems that it should be set equal to 2x(nodes)x(ppn)x(pmem) . 
#PBS -m bea specifies that PBS should send email when the job begins (b), and when it ends (e) or 
aborts (a). 

#PBS - j oe specifies that the output from stdout and stderr will be merged, intermixed, as stdout. 
This example assumes that the executable dscat is located in the same directory where the code is to 
execute and write its output. If ddscat is located in another directory, simply give the full pathname, 
to it. The qsub command is used to submit the PBS job: 

% qsub pbs. submit 

As the calculation proceeds, the usual output files will be written to this directory: for each wave- 
length, target size, and target orientation, there will be a file viaarbbkccc . sea, where aa=00, 01, 02, 
... specifies the wavelength, bb=00, 01, 02, ... specifies the target size, and cc=000, 001, 002, ... spec- 
ifies the orientation. For each wavelength and target size there will also be a file waarbbori . avg 
with orientationally-averaged quantities. Finally, there will also be tables qtable, and qtable2 with 
orientationally-averaged cross sections for each wavelength and target size. 

In addition, each processor employed will write to its own log file ddscat . loqjinn, where nnn=000, 
001, 002, ... . These files contain information concerning cpu time consumed by different parts of the cal- 
culation, convergence to the specified error tolerance, etc. If you are uncertain about how the calculation 
proceeded, examination of these log files is recommended. 

25 Graphics and Postprocessing 

25.1 IDL 

At present, we do not offer a comprehensive package for DDSCAT data postprocessing and graphical 
display in IDL. However, there are several developments worth mentioning: First, we offer several output 
capabilities from within DDSCAT: ASCII (see flaB, FORTRAN unformatted binary (see ^Ta2] i. and 
netCDF portable binary (see iil0.3b . Second, we offer several skeleton IDL utilities: 

• bhmie.pro is our translation to IDL of the popular Bohren-Huffman code which calculates 
efficiencies for spherical particles using Mie theory. 

• readbin.pro reads the FORTRAN unformatted binary file written by routine writebin.f. 
The variables are stored in a common block. 

• readnet . pro reads NetCDF portable binary file and should be the method of choice for IDL 
users. It offers random data access. 

• mie . pro is an example of an interface to binary files and the Bohren-Huffman code. It plots a 
comparison of DDSCAT results with scattering by equivalent radius spheres. 

^ ^http : / / www . openpbs . org| 
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At present the IDL code is experimental. 

26 Miscellanea 

Additional source code, refractive index files, etc., contributed by users will be located in the directory 
DDA/misc. These routines and files should be considered to be not supported by Draine and Flatau - 
caveat receptor! These routines and files should be accompanied by enough information (e.g., comments 
in source code) to explain their use. 

27 Finale 

This User Guide is somewhat inelegant, but we hope that it will prove useful. The structure of the 
ddscat . par file is intended to be simple and suggestive so that, after reading the above notes once, 
the user may not have to refer to them again. 

The file rel_notes in DDA/doc hsts known bugs in DDSCAT 6.1 . Up-to-date release notes will 
be available at the DDSCAT web site, 

http : / / www . astro . princet on . edu/[ ^draine/DDSCAT . html 

Users are encouraged to provide B. T. Draine (draine@astro.princeton.edu) with their 
email address; bug reports, and any new releases of DDSCAT, will be made known to those who do! 

P. J. Flatau maintains the WWW page "SCATTERLIB - Light Scattering Codes Library" with URL 
| http ://at ol.ucsd. edu/^pf latau| The SCATTERLIB Internet site is a library of light scat- 
tering codes. Emphasis is on providing source codes (mostly FORTRAN). However, other information 
related to scattering on spherical and non-spherical particles is collected: an extensive list of references 
to light scattering methods, refractive index, etc. This URL page contains section on the discrete dipole 
approximation. 

Concrete suggestions for improving DDSCAT (and this User Guide) are welcomed. If you wish to 
cite this User Guide, we suggest the following citation: 

Draine, B.T., & Flatau, RJ. 2004, "User Guide for the Discrete Dipole Approximation Code DDSCAT 
(Version 6.1)", |http://arxiv.org/abs/astro-ph/0409262| 

Finally, the authors have one special request: We would very much appreciate preprints and (espe- 
cially!) reprints of any papers which make use of DDSCAT! 
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A Understanding and Modifying ddscat.par 

In order to use DDSCAT to perform the specific calculations of interest to you, it will be necessary to 
modify the ddscat . par file. Here we list the sample ddscat . par file, followed by a discussion of 
how to modify this file as needed. Note that all numerical input data in DDSCAT is read with free-format 
READ (IDEV, *) . . . statements. Therefore you do not need to worry about the precise format in which 
integer or floating point numbers are entered on a line. The crucial thing is that lines in ddscat . par 
containing numerical data have the correct number of data entries, with any informational comments 
appearing after the numerical data on a given line. 

' =================== Parameter file ===================' 

'**** PRELIMINARIES ****' 

'NOTORQ'= CMT0RQ*6 (DOTORQ, NOTORQ) — either do or skip torque calculations 

'PBCGST'= CMDS0L*6 (PBCGST, PETRKP) — select solution method 

'GPFAFT'= CMETHD*6 (GPFAFT, FFTWJ, CONVEX) 

'LATTDR'= CALPHA*6 (LATTDR, SCLDR) 

'NOTBIN'= CBINFLAG (ALLBIN, ORIBIN, NOTBIN) 

'NOTCDF'= CNETFLAG (ALLCDF, ORICDF, NOTCDF) 

' RCTNGL' = CSHAPE*6 (FRMFIL, ELLIPS, CYLNDR, RCTNGL, HEXGON, TETRAH, UNICYL, . . . ) 
32 24 16 = shape parameters PARI, PAR2, PAR3 

1 = NCOMP = number of dielectric materials 

' TABLES' = CDIEL*6 (TABLES, H20ICE, H20LIQ; if TABLES, then filenames follow...) 
'diel.tab' = name of file containing dielectric function 
'**** CONJUGATE GRADIENT DEFINITIONS ****' 

= INIT (TO BEGIN WITH | X0> = 0) 

l.OOe-5 = TOL = MAX ALLOWED (NORM OF | G>=AC I E>-ACA | X> ) / (NORM OF AC|E>) 

Angular resolution for calculation of <cos>, etc. 
0.5 = ETASCA (number of angles is proportional to [ (3+x) /ETASCA] " 2 ) 
'**** Wavelengths (micron) ****' 

6.283185 6.283185 1 ' INV = wavelengths ( first , last , how many , how=LIN, INV, LOG) 
'**** Effective Radii (micron) **** ' 

2 2 1 'LIN' = eff. radii (first, last, how many, how=LIN, INV, LOG) 

Define Incident Polarizations ****' 
(0,0) (l.,0.) (0.,0.) = Polarization state eOl (k along x axis) 
2 = lORTH (=1 to do only pol . state eOl; =2 to also do orth. pol. state) 

1 = IWRKSC (=0 to suppress, =1 to write ".sea" file for each target orient. 
'**** Prescribe Target Rotations ****' 

0. 0. 1 = BETAMI, BETAMX, NBETA (beta=rotation around al) 

0. 90. 3 = THETMI, THETMX, NTHETA (theta=angle between al and k) 

0. 0. 1 = PHIMIN, PHIMAX, NPHI (phi=rotation angle of al around k) 

'**** Specify first IWAV, IRAD, lORI (normally 0) ****' 

= first IWAV, first IRAD, first lORI (0 to begin fresh) 

'**** Select Elements of S_ij Matrix to Print ****' 

6 = NSMELTS = number of elements of S_ij to print (not more than 9) 
11 12 21 22 31 41 = indices ij of elements to print 
'**** Specify Scattered Directions ****' 

0. 0. 180. 10 = phi, thetan_min, thetan_max, dtheta (in degrees) for plane A 
90. 0. 180. 10 = phi, ... for plane B 
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Lines Comments 

1-2 comment lines - no need to change. 

3 NOTORQ if torque calculation is not required; 
DOTORQ if torque calculation is required. 

4 PBCGST is recommended; other option is PETRKP (see ^T4b . 

5 GPFAFT is supplied as default, but FFTW2 1 is recommended if DDSCAT has been compiled with 
FFTW support (see iii j6.2.5l[T3] ); only other option is CONVEX (valid only on a Convex computer). 

6 LATTDR (Draine & Goodman [2000] Lattice Dispersion Relation is only option now supported) 
(see gB. 

7 ALLBIN for full unformatted binary dump ( ^10.21 ); 

ORIBIN for unformatted binary dump of orientational averages only; 
NOTBIN for no unformatted binary output. 

8 ALLCDF for output in netCDF format (must have netCDF option enabled; cf. j)10.3l l; 
ORICDF for orientational averages in netCDF format (must have netCDF option enabled). 
NOTCDF for no output in netCDF format; 

9 specify choice of target shape (see fJT9]for description of options RCTNGL, ELLIPS, TETRAH, ...) 

10 shape parameters SHPARl, SHPAR2, SHPAR3, ... (see 

1 1 number of different dielectric constant tables ( ifT2b . 

12 TABLES - dielectric function to be provided via input table. 

13 name(s) of dielectric constant table(s) (one per line). 

14 comment line - no need to change. 

15 is recommended value of parameter INIT. 

16 TOL = error tolerance h: maximum allowed value of \ A^E — A^ AP\/\A^ E\ [see eq.lfTSlll. 

17 comment line - no need to change. 

18 ETASCA - parameter ?/ controlling angular averages ( ^22^ . 

19 comment line - no need to change. 

20 A - first, last, how many, how chosen. 

21 comment line - no need to change. 

22 Ooff - first, last, how many, how chosen. 

23 comment line - no need to change. 

24 specify x,y,z components of (complex) incident polarization epi (§21]) 

25 lORTH = 1 to do one polarization state only; 

2 to do second (orthogonal) incident polarization as well. 

26 IWRKSC = to suppress writing of ".sea" files; 
2 to enable writing of ".sea" files. 

27 comment line - no need to change. 

28 P (see ifTTl ) - first, last, how many . 

30 Q - first, last, how many. 

31 $ - first, last, how many. 

32 comment line - no need to change. 

33 IWAVO IRADO lORlO - specify Starting values of integers I WAV IRAD I ORI (normally 

34 comment line - no need to change. 

35 A^5 = number of scattering matrix elements (must be < 9) 

36 indices ij of Ns elements of the scattering matrix Sij 

37 comment line - no need to change. 

38 (j>s for first scattering plane, 9s,min, (^s.max, how many Og values; 
39,... (j)s for 2nd,... scattering plane, ... 
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B wxxxr3;_yori . avg Files 

The file wOOOrOOori . avg contains the resuhs for the first wavelength (wOOO) and first target radius 
(rOO) averaged over orientations (ori . avg). The wOOOrOOori . avg file generated by the sample 
calculation should look like the following: 

DDSCAT DDSCAT 6.1 [05.10.18] 

TARGET Rectangular prism; NX,NY,NZ= 32 24 15 

GPFAFT method of solution 

LATTDR prescription for polarizabilies 

RCTNGL shape 

12288 = NATO = number of dipoles 

AEFF= 2.00000 = effective radius (physical units) 
WAVE= 6.28319 = wavelength (physical units) 
K*AEFF= 2.00000 = 2*pi*aeff /lambda 

n= ( 1.3300 , 0.0100), eps.= ( 1.7688 , 0.0266) | m | kd= 0.1858 for subs. 1 
TOL= l.OOOE-05 = error tolerance for CCG method 
NAVG= 603 = (theta,phi) values used in comp . of Qsca, g 
( 1.00000 0.00000 0.00000) = target axis Al in Target Frame 
( 0.00000 1.00000 0.00000) = target axis A2 in Target Frame 
( 0.13971 0.00000 0.00000) = k vector (latt. units) in Lab Frame 
( 0.00000, 0.00000) ( 1.00000, 0.00000) ( 0.00000, O.00000)=inc.pol.vec. 1 in LF 
( 0.00000, 0.00000) ( 0.00000, 0.00000) ( 1.00000, . 00000) =inc .pol .vec . 2 in LF 
0.000 0.000 = beta_min, beta_max ; NBETA = 1 
0.000 90.000 = theta_min, theta_max; NTHETA= 3 
0.000 0.000 = phi_min, phi_max ; NPHI = 1 
0.5000 = ETASCA = param. controlling # of scatt. dirs used to calculate <cos> etc. 
Results averaged over 3 target orientations 
and 2 incident polarizations 
Qext Qabs Qsca g(l)=<cos> <cos"2> Qbk Qpha 

J0=1: 8.6315E-01 8.1919E-02 7.8123E-01 6.5936E-01 5.5097E-01 8.6249E-03 9.5916E-01 
J0=2: 5.9757E-01 6.5818E-02 5.3175E-01 7.0839E-01 5.7267E-01 6.6347E-03 8.9156E-01 
mean: 7.3036E-01 7.3868E-02 6.5649E-01 6.7921E-01 5.7267E-01 7.6298E-03 9.2536E-01 
Qpol= 2.6558E-01 dQpha= 6.7601E-02 

Qsca*g(l) Qsca*g(2) Qsca*g(3) iter mxiter Nsca 



J0= 


= 1 : 


5 


1511E 


-01 2 . 


0865E-02 - 


2 . 


3925E 


-08 




5 


10000 


603 














J0= 


=2 : 


3 


7669E 


-01 3 . 


3540E-02 - 


2 . 


2460E 


-08 




4 


10000 


603 














mean : 


4 


4590E 


-01 2 . 


7202E-02 - 


2 . 


3192E 


-08 
































■k -k 


Mueller matrix elements 


:or 


selected scattering directions 










theta 


phi 




Pol . 




3_11 




S_12 






S_21 




C 


_2 2 




S_31 




c 


_41 







.0 








. 


11130 


3 


. 982E+00 


4 


.433E 


-01 


4 


. 433E- 


01 


3 . 


982E+00 


-2 


. 173E- 


08 


-2 . 


225E- 


08 


10 


.0 








. 


09795 


3 


. 774E+00 


3 


. 697E 


-01 


3 


. 697E- 


01 


3 . 


774E+00 


-9 


. 126E- 


10 


-2 . 


229E- 


08 


20 


.0 








. 


06502 


3 


. 211E + 00 


2 


.088E 


-01 


2 


. 088E- 


01 


3. 


211E+00 


9 


. 185E- 


09 


-1 . 


218E- 


09 


30 


.0 








. 


01089 


2 


. 475E+00 


2 


. 695E 


-02 


2 


. 695E- 


02 


2 . 


475E+00 


3 


. 537E- 


09 


-4 . 


172E- 


08 


40 


.0 








. 


06563 


1 


. 745E+00 


-1 


. 145E 


-01 


-1 


. 145E- 


01 


1 . 


745E+00 


9 


. 433E- 


09 


3 . 


165E- 


09 


50 


.0 








. 


16356 


1 


. 138E+00 


-1 


.861E 


-01 


-1 


. 861E- 


01 


1 . 


138E+00 


-3 


. 586E- 


08 


9. 


124E- 


09 


60 


.0 








. 


27608 


6 


. 918E-01 


-1 


. 910E 


-01 


-1 


. 910E- 


01 


6. 


918E-01 


2 


. 853E- 


08 


2 . 


947E- 


09 


70 


.0 








. 


38508 


3 


. 932E-01 


-1 


. 514E 


-01 


-1 


. 514E- 


01 


3 . 


932E-01 


7 


. 071E- 


08 


8 . 


935E- 


09 


80 


.0 








. 


45639 


2 


. 066E-01 


-9 


. 431E 


-02 


-9 


. 431E- 


02 


2 . 


066E-01 


-1 


.103E- 


08 


-2 . 


167E- 


08 


90 


.0 








. 


43728 


9 


. 729E-02 


-4 


.254E 


-02 


-4 


. 254E- 


02 


9. 


729E-02 


1 


.201E- 


08 


-1 . 


803E- 


08 


100 


.0 








. 


24465 


4 


. 032E-02 


-9 


. 864E 


-03 


-9 


. 864E- 


03 


4 . 


032E-02 


1 


. 691E- 


08 


5 . 


428E- 


09 


110 


.0 








. 


09357 


2 


. 051E-02 


1 


. 919E 


-03 


1 


. 919E- 


03 


2 . 


051E-02 


-6 


. 085E- 


10 


2 . 


074E- 


09 


120 


.0 








. 


03725 


2 


. 702E-02 


1 


.006E 


-03 


1 


. 006E- 


03 


2 . 


702E-02 


-3 


. 231E- 


09 


3 . 


086E- 


09 


130 


.0 








. 


03659 


4 


. 849E-02 


-1 


. 774E 


-03 


-1 


. 774E- 


03 


4 . 


849E-02 


-6 


. 712E- 


10 


5 . 


339E- 


10 


140 


.0 








. 


00643 


7 


. 271E-02 


-4 


. 672E 


-04 


-4 


. 672E- 


04 


7 . 


271E-02 


-6 


. 046E- 


09 


1 . 


116E- 


10 


150 


.0 








. 


04894 


9 


. 069E-02 


4 


.439E 


-03 


4 


. 439E- 


03 


9. 


069E-02 


2 


. 738E- 


09 


1 . 


609E- 


09 


160 


.0 








. 


09716 


9 


. 955E-02 


9 


. 672E 


-03 


9 


. 672E- 


03 


9. 


955E-02 


4 


. 978E- 


09 


-9. 


134E- 


10 


170 


.0 








. 


12616 


1 


. 007E-01 


1 


.271E 


-02 


1 


. 271E- 


02 


1 . 


007E-01 


-5 


. 866E- 


10 


-2 . 


023E- 


10 


180 


.0 








. 


13043 


9 


. 588E-02 


1 


.251E 


-02 


1 


. 251E- 


02 


9. 


588E-02 


-2 


.298E- 


10 


3 . 


624E- 


10 





.0 


90 





. 


11130 


3 


. 982E+00 


-4 


.433E 


-01 


-4 


. 433E- 


01 


3. 


982E+00 


-1 


. 986E- 


08 


-1 . 


817E- 


08 
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10 





90 








12306 


3 


872E+00 


-4 


764E- 


01 


-4 


763E- 


01 


3 


872E+00 


-1 


299E- 


02 


-1 


096E- 


02 


20 





90 








15885 


3 


553E+00 


-5 


642E- 


01 


-5 


640E- 


01 


3 


553E+00 


-2 


416E- 


02 


-2 


003E- 


02 


30 





90 








22005 


3 


064E+00 


-6 


738E- 


01 


-6 


734E- 


01 


3 


063E+00 


-3 


172E- 


02 


-2 


559E- 


02 


40 





90 








30823 


2 


466E+00 


-7 


598E- 


01 


-7 


593E- 


01 


2 


465E+00 


-3 


428E- 


02 


-2 


658E- 


02 


50 





90 








42360 


1 


842E+00 


-7 


799E- 


01 


-7 


794E- 


01 


1 


841E+00 


-3 


144E- 


02 


-2 


298E- 


02 


60 





90 








56195 


1 


270E+00 


-7 


135E- 


01 


-7 


132E- 


01 


1 


269E+00 


-2 


421E- 


02 


-1 


598E- 


02 


70 





90 








71056 


8 


064E-01 


-5 


729E- 


01 


-5 


728E- 


01 


8 


060E-01 


-1 


487E- 


02 


-7 


702E- 


03 


80 





90 








84571 


4 


703E-01 


-3 


976E- 


01 


-3 


977E- 


01 


4 


701E-01 


-6 


076E- 


03 


-3 


938E- 


04 


90 





90 








93824 


2 


500E-01 


-2 


345E- 


01 


-2 


345E- 


01 


2 


498E-01 


1 


803E- 


04 


4 


398E- 


03 


100 





90 








96927 


1 


190E-01 


-1 


155E- 


01 


-1 


153E- 


01 


1 


186E-01 


3 


172E- 


03 


6 


275E- 


03 


110 





90 








93387 


5 


089E-02 


-4 


795E- 


02 


-4 


740E- 


02 


5 


023E-02 


3 


379E- 


03 


5 


851E- 


03 


120 





90 








76425 


2 


524E-02 


-2 


OOOE- 


02 


-1 


919E- 


02 


2 


436E-02 


1 


957E- 


03 


4 


220E- 


03 


130 





90 








47050 


2 


657E-02 


-1 


339E- 


02 


-1 


250E- 


02 


2 


563E-02 


1 


204E- 


04 


2 


395E- 


03 


140 





90 








31152 


4 


231E-02 


-1 


390E- 


02 


-1 


312E- 


02 


4 


150E-02 


-1 


280E- 


03 


9 


816E- 


04 


150 





90 








22721 


6 


247E-02 


-1 


461E- 


02 


-1 


407E- 


02 


6 


192E-02 


-1 


882E- 


03 


1 


594E- 


04 


160 





90 








17286 


8 


023E-02 


-1 


404E- 


02 


-1 


376E- 


02 


7 


995E-02 


-1 


712E- 


03 


-1 


589E- 


04 


170 





90 








14100 


9 


187E-02 


-1 


299E- 


02 


-1 


292E- 


02 


9 


180E-02 


-9 


922E- 


04 


-1 


540E- 


04 


180 





90 








13043 


9 


588E-02 


-1 


251E- 


02 


-1 


251E- 


02 


9 


588E-02 


2 


708E- 


09 


9 


289E- 


10 
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C ^fjxxxryykzzz .sea Files 

The wOOOrOOkOOO.sca file contains the resuhs for the first wavelength (wOOO), first target radius 
(rOO), and first orientation (kOOO). The wOOOrOOkOOO . sea file created by the sample calculation 
should look like the following; 

DDSCAT DDSCAT 6.1 [05.10.18] 

TARGET Rectangular prism; NX,NY,NZ= 32 24 16 

GPFAFT method of solution 

LATTDR prescription for polarizabilies 

RCTNGL shape 

12288 = NATO = number of dipoles 

AEFF= 2.00000 = effective radius (physical units) 
WAVE= 6.28319 = wavelength (physical units) 
K*AEFF= 2.00000 = 2*pi*aeff /lambda 

n= ( 1.3300 , 0.0100), eps.= ( 1.7688 , 0.0266) I m | kd= 0.1858 for subs. 1 
TOL= l.OOOE-05 = error tolerance for COG method 
NAVG= 603 = (theta,phi) values used in comp. of Qsca,g 
( 1.00000 0.00000 0.00000) = target axis Al in Target Frame 
( 0.00000 1.00000 0.00000) = target axis A2 in Target Frame 
( 0.13971 0.00000 0.00000) = k vector (latt. units) in TF 

( 0.00000, 0.00000) ( 1.00000, 0.00000) ( 0.00000, O.00000)=inc.pol.vec. 1 in TF 
( 0.00000, 0.00000) ( 0.00000, 0.00000) ( 1.00000, O.00000)=inc.pol.vec. 2 in TF 
BETA = 0.000 = rotation of target around Al 
THETA= 0.000 = angle between Al and k 

PHI = 0.000 = rotation of Al around k 
0.5000 = ETASCA = param. controlling # of scatt. dirs used to calculate <cos> etc. 

Qext Qabs Qsca g(l)=<cos> <cos"2> Qbk Qpha 

J0=1: 9.0975E-01 8.6073E-02 8.2368E-01 6.4091E-01 5.8433E-01 2.5118E-02 9.7721E-01 
J0=2 : 6.9871E-01 7.0481E-02 6.2823E-01 6.5890E-01 6.0513E-01 1.9109E-02 9.1273E-01 
mean: 8.0423E-01 7.8277E-02 7.2595E-01 6.4870E-01 5.9333E-01 2.2114E-02 9.4497E-01 
Qpol= 2.1104E-01 dQpha= 6.4479E-02 

Qsca*g(l) Qsca*g(2) Qsca*g(3) iter mxiter Nsca 



J0= 


= 1 : 


5 


2791E 


-01 2 . 


1428E-08 - 


4 . 


2616E 


-08 




5 


10000 


603 














J0= 


= 2 : 


4 


1394E 


-01 -5. 


8902E-08 


2 . 


9659E 


-08 




4 


10000 


603 














mean : 


4 


7092E 


-01 -1. 


8737E-08 - 


6. 


4789E 


-09 
































* * 


Mueller 


Tiatrix elements 


for 


selected scattering directions 










theta 


phi 




Pol . 




S_ll 




S_12 






S_21 






S_2 2 


c 


_31 




c 


_41 
















. 


09765 


4 


. 234E+00 


4 


. 134E 


-01 


4 


. 134E- 


01 


4 


. 234E + 00 


-2 . 


913E- 


09 


-1 . 


334E- 


09 


10 











. 


08807 


4 


. 052E+00 


3 


.569E 


-01 


3 


. 569E- 


01 


4 


. 052E+00 


-2 . 


714E- 


08 


3 . 


066E- 


08 


20 











. 


05875 


3 


. 544E+00 


2 


. 082E 


-01 


2 


. 082E- 


01 


3 


. 544E+00 


-1 . 


112E- 


08 


2 . 


186E- 


08 


30 











. 


00802 


2 


. 818E+00 


2 


. 261E 


-02 


2 


. 261E- 


02 


2 


. 818E+00 


-4 . 


716E- 


08 


-1 . 


355E- 


07 


40 











. 


06640 


2 


. 023E+00 


-1 


.343E 


-01 


-1 


. 343E- 


01 


2 


. 023E+00 


-1 . 


463E- 


09 


-6. 


662E- 


08 


50 











. 


16610 


1 


. 301E+00 


-2 


. 160E 


-01 


-2 


. 160E- 


01 


1 


. 301E + 00 


4 . 


015E- 


08 


1 . 


013E- 


07 


60 











. 


28816 


7 


. 438E-01 


-2 


. 143E 


-01 


-2 


. 143E- 


01 


7 


. 438E-01 


3. 


558E- 


08 


-7 . 


154E- 


09 


70 











. 


41587 


3 


. 746E-01 


-1 


.558E 


-01 


-1 


. 558E- 


01 


3 


. 746E-01 


2 . 


146E- 


09 


4 . 


203E- 


09 


80 











. 


49780 


1 


. 624E-01 


-8 


. 084E 


-02 


-8 


. 084E- 


02 


1 


. 624E-01 


-3. 


713E- 


08 


-6. 


OOOE- 


09 


90 











. 


39313 


5 


. 674E-02 


-2 


. 231E 


-02 


-2 


. 231E- 


02 


5 


. 674E-02 


-1 . 


134E- 


09 


6. 


122E- 


09 


100 











. 


34236 


1 


. 474E-02 


5 


. 046E 


-03 


5 


. 046E- 


03 


1 


. 474E-02 


-8 . 


243E- 


09 


4 . 


109E- 


09 


110. 











. 


39371 


1 


. 130E-02 


4 


. 449E 


-03 


4 


. 449E- 


03 


1 


. 130E-02 


-4 . 


827E- 


10 


-3 . 


227E- 


09 


120 











. 


25972 


3 


. 555E-02 


-9 


. 234E 


-03 


-9 


. 234E- 


03 


3 


.555E-02 


-4 . 


603E- 


09 


-9. 


526E- 


09 


130 











. 


23736 


8 


. 116E-02 


-1 


. 926E 


-02 


-1 


. 926E- 


02 


8 


. 116E-02 


-1 . 


637E- 


09 


3 . 


782E- 


09 


140 











. 


11691 


1 


. 388E-01 


-1 


. 622E 


-02 


-1 


. 622E- 


02 


1 


.388E-01 


1 . 


650E- 


10 


-2 . 


609E- 


09 


150 











. 


00665 


1 


. 957E-01 


-1 


.301E 


-03 


-1 


. 301E- 


03 


1 


. 957E-01 


-5. 


184E- 


09 


-2 . 


104E- 


09 


160 











. 


07318 


2 


. 409E-01 


1 


.763E 


-02 


1 


. 763E- 


02 


2 


. 409E-01 


3. 


098E- 


09 


-3 . 


475E- 


09 


170 











. 


12034 


2 


. 687E-01 


3 


.233E 


-02 


3 


.233E- 


02 


2 


. 687E-01 


-4 . 


327E- 


09 


4 . 


008E- 


10 


180 











. 


13587 


2 


. 779E-01 


3 


. 776E 


-02 


3 


. 776E- 


02 


2 


. 779E-01 


-3. 


491E- 


10 


-1 . 


861E- 


10 








90 





. 


09765 


4 


. 234E+00 


-4 


. 134E 


-01 


-4 


. 134E- 


01 


4 


. 234E + 00 


-3. 


785E- 


08 


-1 . 


094E- 


09 


10 





90 





. 


10919 


4 


. 116E+00 


-4 


. 494E 


-01 


-4 


. 494E- 


01 


4 


. 116E+00 


-4 . 


023E- 


08 


-1 . 


834E- 


09 
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20 





90 








14431 


3 


770E+00 


-5 


440E- 


01 


-5 


440E- 


01 


3 


770E+00 


-2 


545E- 


08 


6 


757E- 


09 


30 





90 








20431 


3 


226E+00 


-6 


591E- 


01 


-6 


591E- 


01 


3 


226E+00 


-1 


610E- 


08 


6 


382E- 


09 


40 





90 








29061 


2 


548E+00 


-7 


405E- 


01 


-7 


405E- 


01 


2 


548E+00 


-2 


430E- 


08 


-6 


944E- 


09 


50 





90 








40308 


1 


829E+00 


-7 


373E- 


01 


-7 


373E- 


01 


1 


829E+00 


-1 


432E- 


08 


1 


029E- 


08 


60 





90 








53705 


1 


171E+00 


-6 


291E- 


01 


-6 


291E- 


01 


1 


171E+00 


-1 


211E- 


08 


1 


836E- 


08 


70 





90 








67862 


6 


517E-01 


-4 


422E- 


01 


-4 


422E- 


01 


6 


517E-01 


-1 


765E- 


08 


3 


355E- 


10 


80 





90 








80009 


2 


991E-01 


-2 


393E- 


01 


-2 


393E- 


01 


2 


991E-01 


-3 


285E- 


09 


7 


309E- 


09 


90 





90 








85380 


9 


905E-02 


-8 


457E- 


02 


-8 


457E- 


02 


9 


905E-02 


-2 


598E- 


09 


6 


533E- 


09 


100 





90 








60725 


1 


530E-02 


-9 


294E- 


03 


-9 


294E- 


03 


1 


530E-02 


-4 


388E- 


10 


-6 


103E- 


10 


110 





90 








20216 


1 


053E-02 


-2 


128E- 


03 


-2 


128E- 


03 


1 


053E-02 


-1 


217E- 


10 


-6 


628E- 


10 


120 





90 








51167 


5 


332E-02 


-2 


728E- 


02 


-2 


728E- 


02 


5 


332E-02 


-3 


996E- 


11 


-1 


717E- 


09 


130 





90 








44322 


1 


171E-01 


-5 


188E- 


02 


-5 


188E- 


02 


1 


171E-01 


1 


269E- 


09 


-1 


079E- 


09 


140 





90 








34094 


1 


800E-01 


-6 


136E- 


02 


-6 


136E- 


02 


1 


800E-01 


1 


652E- 


09 


-3 


556E- 


09 


150 





90 








25164 


2 


285E-01 


-5 


749E- 


02 


-5 


749E- 


02 


2 


285E-01 


1 


863E- 


09 


3 


855E- 


10 


160 





90 








18685 


2 


587E-01 


-4 


833E- 


02 


-4 


833E- 


02 


2 


587E-01 


2 


355E- 


09 


6 


633E- 


12 


170 





90 








14849 


2 


736E-01 


-4 


063E- 


02 


-4 


063E- 


02 


2 


736E-01 


3 


235E- 


09 


-2 


432E- 


10 


180 





90 








13587 


2 


779E-01 


-3 


776E- 


02 


-3 


776E- 


02 


2 


779E-01 


3 


490E- 


09 


-1 


259E- 


10 



Index 



Qcxt, 9 
Qsca,^ 9 

S'j - scattering amplitude matrix, 42 

Stj - 4x4 Mueller scattering matrix, 42, 43, 45 

S'j - scattering amplitude matrix, 43 

$ - target orientation angle, 19, 30 

G - target orientation angle, 19, 30 

P - target orientation angle, 19, 30 

7] - parameter for choosing scattering angles, 40 

J] - parameter for selection of scattering angles, 40 

ai, 3.2 target vectors, 19 

(pg - scattering angle, 38, 40 

9s - scattering angle, 38, 40 

Ooff, 6, 19 

fij - amplitude scattering matrix, 42 
fij - scattering amplitude matrix, 43 
Qr - radiative torque efficiency vector, 27 
Qpr - radiative force efficiency vector, 27 

absorption efficiency factor Qabs> 9 
averages over scattering angles, 38 

CALLTARGET, 37 
compiling and linking, 13, 17 
conjugate gradient algorithm, 26 
cxfftw_fake.f, 15 

ddscat.log.OOO - output file, 20 

ddscat.out, 20 

ddscat.par 

ALLBIN, 21 

CSHAPE, 33 

diel.tab, 19 

DOTORQ, 27 

ETASCA, 28 

FFTW21, 15 

GKDLDR, 22 

GPFAFT, 19 

lORTH, 42, 46 

IWRKSC, 21 

LATTDR, 19, 22 

NOTBIN, 19 

NOTCDF, 19 

ORIBIN, 21 

PBCGST, 27 

PETRKP, 27 

SHPAR1,SHPAR2,..., 33 

TOL, 27 

ddscat.par - BETAMI,BETAMX,NBETA, 31 
ddscat.par - PHIMIN,PHIMAX,NPHI, 31 



ddscat.par - THETMI,THETMX,NTHETA, 31 

ddscat.par parameter file, 19 

diel.tab - see ddscat.par, 19 

dielectric function of target material, 23 

dielectric functions, multiple, 23 

dielectric medium, 9 

differential scattering cross section dCsca/dH., 45 
DOTORQ, 27 

effective radius aofr, 6, 19 
ELLIPS, 33 
error tolerance, 27 
ESELE 49 
ETASCA, 28 

extinction efficiency factor Qoxt, 9 

FFLAGSslamcl flag in Makefile, 13, 15 
FFT algoriflim, 24 
FFTW, 14, 15, 24, 49 
FFTW21 - see ddscat.par, 15 
fortran compiler, 14 

optimization, 13, 15 

GKDLDR, 22 
GPFA, 15, 24 

GPFAFT - see ddscat.par, 19 
GPFAPACK, 49 

IDL,21,48 

bhmie.pro, 48 

mie.pro, 48 

readbin.pro, 48 

readnet.pro, 48 
IDVERR, 14 
IDVOUT, 14 
lORTH, 46 
lORTH parameter, 42 
iterative algorithm, 26 
IWRKSC - see ddscat.par, 21 

Lab Frame (LF), 29 
LATTDR, 22 

LATTDR - see ddscat.par, 19 

Makefile, 14, 17 

memory requirements, 28 

MPI - Message Passing Interface, 14, 46 

code execution, 47 

Makefile, 47 

mpif.h, 47 

mpisubs.f, 47 
mtable - output file, 20 



58 



INDEX 



59 



Mueller matrix for scattering, 42, 43, 45 
MXCOMP, 28 
MXNX,MXNY,MXNZ, 28 

netCDF, 14,21 
netCDF capability 

disabled, 16 

enabled, 17 
NOTBIN - see ddscat.par, 19 
NOTCDF - see ddscat.par, 19 

orientational averaging, 3 1 

nonrandomly-oriented targets, 32 
randomly-oriented targets, 32 

orientational sampling in (3, G, and $,31 

PBCGST algorithm, 27 

PETRKP algorithm, 27 

phase lag efficiency factor Qpha, 9 

PIM package, 27 

polarizabilities 

GKDLDR, 22 

LATTDR, 22 
polarization - elliptical, 42 
polarization of incident radiation, 38 
polarization of scattered radiation, 44, 45 

qtable - output file, 20 
qtable2 - output file, 20 
qtable2 file, 20 

radiative force efficiency vector Qpr, 27 
radiative torque efficiency vector Qr, 27 
REFICE, 49 

refractive index of target material, 23 
refractive indices, multiple, 23 
REFWAT, 49 

relative dielectric function, 9 
relative refractive index, 10 

scattering - angular averages, 38 
scattering angles Bg, (f>s, 40 
scattering by tilted cube, 38 
scattering directions, 38 
scattering efficiency factor Qsca, 9 
SCATTERLIB, 49 
size parameter x = kacs, 6 
SLAMCl - do not optimize, 13, 15 
source code (downloading of), 1 1 
Stokes vector (/, Q, U, V), 43 

Target Frame (TF), 30 
target generation, 32 
target orientation, 29, 30 
target routines: modifying, 37 



target shape options 

ANIELL, 33 

ANIREC, 33 

BLOCKS, 35 

CONELL, 33 

CYLCAP 33 

CYLNDR, 33 

DW1996, 36 

ELLIPS, 33 

FRMFIL, 37 

HEXGON, 34 

LYRSLB, 37 

NSPHER, 36 

PR1SM3, 36 

RCTNGL, 33 

TETRAH, 34 

THRAEL, 35 

THRELL, 35 

TWOAEL, 35 

TWOELL, 35 

TWOSPH, 34 

UNICYL, 33 
target. out - output file, 38 
TIMEIT subroutine, 14 

wOOOrOOkOOO.sca- output file, 20, 21 
wOOOrOOori.avg - output file, 20 
wOOOrOOori.avg file - output file, 20 
winzip, 12 
WRITENET, 17 
wri tenet. f, 17 
writenet_fake.f, 16 



